Shareware Grab Bag/Shareware Grab Bag.iso/090/byt0486a.arc/ATOMTLBX.LBR/MANUAL.DOC

Text File | 1986-04-11 | 122KB | 3,444 lines

                         ATOMCC USER MANUAL
                         ------ ---- ------

                    For Micro-Computers on MSDOS

                            Version 7.10

                                 by


                            Y. F. Chang























Copyright (C) 1983, 1985 Y. F. Chang

MSDOS is a registered trademark of Microsoft Corp.


                              - i -






                    Table of Contents


Chapter 1 Introduction                                        1

   1.1 The Major Advancements in this Version                 1
   1.2 Purpose and Requirements of the Translator             1
   1.3 Applicability                                          2
   1.4 System Overview                                        4

      1.4.1 The Translator, ATOMCC                            4
      1.4.2 Tha Object Program, ATSPGM                        4

   1.5 Purpose of the User Manual                             5
   1.6 Acknowledgements                                       6
   1.7 References                                             6

Chapter 2 For New Users                                       8

   2.1 Task of the ATOMCC Translator                          8
   2.2 Using the ATOMCC system                                9

      2.2.1 Step 1 - edit ODEINP                              10
      2.2.2 Step 2 - Run ATOMCC                               14
      2.2.3 Step 3 and 4 - Compile and link ATSPGM            15
      2.2.4 Step 5 - Prepare the data                         15
      2.2.5 Step 6 - Run ATSPGM                               16

   2.3 Output at equally spaced points                        17

      2.3.1 ZEROT - Stopping at roots of variables            17
      2.3.2 MSTIFF=20,21,22 - Stiff problems.                 18

Chapter 3 How to Use ----                                     20

   3.1 Solving your problem                                   20

      3.1.1 Translator file, ODEINP                           20
      3.1.2 Translator file, the terminal                     22
      3.1.3 Translator file, ATSPGM                           23
      3.1.4 DATA input file                                   27
      3.1.5 Solution file                                     28
      3.1.6 User files                                        29

   3.2 Using block 1                                          29

      3.2.1 Format for the system of equations                30
      3.2.2 Parameters in the equations                       31
      3.2.3 COPTION DOUBLE - Double-precision ATSPGM          31
      3.2.4 COPTION COMPLX - Complex ATSPGM                   32
      3.2.5 COPTION DOUBLE, COMPLX - Double-complex ATSPGM    33
      3.2.6 COPTION LENVAR=n - Series length                  35
      3.2.7 COPTION DUMP=n - Diagnostic messages              36


                              - ii -



   3.3 Using block 2                                          36

      3.3.1 Subroutine form of ATSPGM                         37
      3.3.2 User declarations                                 37
      3.3.3 Common blocks for user                            38

   3.4 Using block 3                                          38

      3.4.1 Initial conditions                                39
      3.4.2 Parameters in the differential equations          39
      3.4.3 Solve a problem repeatedly                        40
      3.4.4 START, END - Interval of integration              40
      3.4.5 NSTEPS - Number of integration steps, default=40  40
      3.4.6 H - Initial trial stepsize, default = 1.0         41
      3.4.7 ERRLIM - Preset accuracy of the solution          41
      3.4.8 ADJSTF - Error control for stiff problems         41
      3.4.9 LENSER - Length of series used, default=30        42
      3.4.10 MPRINT - Amount of print produced, default=4     42
      3.4.11 DLTXPT - Print point increments, default = 0.0   43
      3.4.12 KTRDCV - Dynamic suppression of CALL RDCV        43
      3.4.13 KPTS - Number of points on complex path          44
      3.4.14 POINTS - Complex path of integration             44
      3.4.15 MSTIFF=10 - Solutions which are entire           44
      3.4.16 MSTIFF=20,21,22 - Stiff problems.                44
         3.4.16.1 Steady-State Stiff Problems                 45

   3.5 Using block 4                                          46

      3.5.1 Automatic printing of output points               46
      3.5.2 User controlled printing of output points         47
      3.5.3 Logarithmic spacing of output points              48
      3.5.4 ZEROT - Stopping & printing at roots of variables 48
      3.5.5 Finding singularities in real solutions           49
      3.5.6 Stopping short of a singularity                   51

   3.6 Editing of ATSPGM                                      51

      3.6.1 TERM - Fast generation of print at output points  51

   3.7 Large systems                                          54
   3.8 Solving ODE's in the complex domain                    54


                               - 1 -








                             Chapter 1

                            Introduction



    This chapter is written to help you become familiar with the
purpose and requirements of the ATOMCC system and with the
organization of this manual.




1.1 The Major Advancements in this Version


    The present version of ATOMCC, 7.10 for micros, contains a
major advancement.  Now, the ATOMCC system will solve stiff
problems.  This represents a significant departure from the central
premise of the ATOMCC system, which is precise error control.  For
non-stiff problems, the user still have the most accurately
controlled numerical method ever developed.  For many problems, the
precision is so good that there is ALMOST global error control.

    For stiff problems, due to the nature of the "approximating"
solution, there cannot be true error control.  Therefore, the
controlling parameter for errors in stiff problems is called
ADJSTF.  It is only meant to be an adjusting constant that can be
loosely refered to as an error control.

    There is also another particularly useful feature in the
current version.  All the dependent variables are now placed into a
temporary two-dimensional array (TMPS) by an EQUIVALENCE
statement.  This allows the user to reference each variable by an
index value.  For a system with x, y, and z as functions of t, the
term y(5) can be also refered to as TMPS(5,2).  Similarly, z(23) =
TMPS(23,3).




1.2 Purpose and Requirements of the Translator


    The ATOMCC system is a tool to be used in the solution of all
initial value problems in ordinary differential equations, (stiff
as well as non-stiff).  It is simple enough to be used by students,
practical enough to be used by engineers, and versatile enough to
be used by research mathematicians.

    The ATOMCC package is delivered on floppy disks.  It uses
Microsoft-FORTRAN77 (a registered trademark of Microsoft Corp.)


                               - 2 -


and works on a micro-computer operating under MS-DOS.  With the
ATOMCC system, you now have in your possession a research tool
whose capabilities far exceed those of standard numerical
integration methods available on main-frame computers.  To be able
to run ATOMCC on your MSDOS micro-computer, you must have the
following hardware and software:-

  -  an MSDOS computer, with an 8087 co-processor;

  -  at least 256K of RAM memory;

  -  two floppy disc drives, or a hard disc drive;

  -  the Microsoft-FORTRAN77 version 3.30.

    A complete system includes the following disc files.

  -  The ATOMCC system files are:-

     ATOMCC.EXE     This is the ATOMCC compiler that reads
                    statements of differential equations and
                    generates an object FORTRAN program called
                    ATSPGM.  The name ATSPGM for the object program
                    file is fixed, but you may change it after it
                    has been written by ATOMCC.

     RDCV.OBJ       This is the ATOMCC subroutine library in
                    single-precision.

     DRDCV.OBJ      This is the ATOMCC subroutine library in
                    double-precision.

     CRDCV.OBJ      This is the ATOMCC subroutine library in
                    complex.

     CDRDCV.OBJ     This is the ATOMCC subroutine library in
                    complex-double.

  -  The Microsoft-FORTRAN77 files are descripted in the Microsoft
     manual.  The relevant files are:- FOR1.EXE, PAS2.EXE,
     LINK.EXE, and FORTRAN.LIB.

    Throughout the discussions in this User Manual, we shall assume
that all of the ATOMCC system files are on the A: floppy-disk
drive, and the Microsoft-FORTRAN77 files are on the B: drive.




1.3 Applicability


  -  The ATOMCC method can solve:

       *  systems of stiff and non-stiff systems of initial value
          problems in ordinary differential equations in which


                               - 3 -


       *  the highest order derivative of each dependent variable
          is given explicitly on the left hand side of an equation,
          whose right hand side has a finite sequence of +, -, *,
          /, **, EXP, SIN, COS, TAN, SINH, COSH, TANH, ALOG, ACOS,
          ASIN, ATAN, or any function which is the solution to a
          differential equation.

  -  The known limitations of ATOMCC are:

       *  the derivatives may be of order at most 6, and

       *  there are at most 900 equations in the system.

  -  ATOMCC can also solve (with manual intervention):

       *  solutions which are polynomials,

       *  singular problems which require the application of
          l'Hopital's rule, or

       *  problems which have catastrophic subtractive errors in
          series generation.

  -  ATOMCC is most attractive for:-

       *  problems with stringent accuracy requirements,

       *  stiff problems,

       *  problems which must be solved repeatedly (such as
          parameter identification), or

       *  quick and easy problems (students' assignments).

  -  The very high order and precise error control used by ATOMCC
     have enabled it to solve many problems which other methods
     were unable to solve.

  -  The ATOMCC compiler allows for the solution of ODE's in the
     complex domain.  This unique capability can be used to explore
     the structure of the singularities in the complex domain of
     non-linear problems.  The analytic information about the
     location and order of singularities in the solution provides
     insight into the behavior of the system.  This method has been
     used to map the first mathematical natural boundary discovered
     in the solution of a nonlinear dynamics problem (7).

  -  The complexity and execution time of ATSPGM depend on the
     number of functions and on the number of multiplications in
     the ODE system, not on the number of equations in the ODE
     system nor on the order of the derivatives involved.  There is
     no penalty for high-order derivatives.

  -  As with all numerical methods, there is no substitute for
     insight into the structure of the ODE system and for the
     application of clever transformations.


                               - 4 -


    Solutions which are entire (have no singularities in the finite
plane), should not be solved using the ATOMCC system.  It is a
total waste of computing power to solve linear problems using
ATOMCC.  This is particularly true for linear 'stiff' problems.

    It can be EASILY show that ALL solutions that are entire can be
solved in quasi-closed forms.  This INCLUDES two-point boundary
value problems!

    Some special circumstances, which rarely occur, are identified
either by ATOMCC or by the series analysis software, and an
appropriate message is produced.  In such cases, the user should
examine the series (using MPRINT=10 in the third block), and seek
the advice of the authors.




1.4 System Overview


1.4.1 The Translator, ATOMCC

    The ATOMCC translator is an ODE-solving compiler written in
FORTRAN.  The ODE system to be solved is written into the ODEINP
input file using conventions discussed in Chapters 2 and 3. The
name ODEINP for the input file is fixed within ATOMCC; you must use
this name.  ATOMCC reads ODEINP and produces a FORTRAN object
program, called ATSPGM.  The name ATSPGM for the object program is
also fixed; you must have a file by this name on your disc even if
it is an empty file.  The numerical solution to the ODE system is
obtained by compiling and executing ATSPGM together with the
library subroutine RDCV.OBJ or one of its variants.

    ATOMCC accepts four blocks of data from ODEINP in which the
user specifies the differential equations, the integration
interval, initial conditions, and various other parameters to be
used in the solution.  The first block is used to specify the
differential equations and commands to ATOMCC.  The second through
fourth blocks are used to insert statements directly into ATSPGM.
The third block is required to specify the integration interval and
the initial conditions.  Detailed guidelines for the use of each
block appear in Chapter 3.


1.4.2 Tha Object Program, ATSPGM

    The ATSPGM object program implements the Taylor series
algorithm for solving initial value problems in ordinary
differential equations.  This Taylor series algorithm is outlined
below.

  -  Initialize method control parameters which may be modified.

  -  Assign initial conditions, the integration interval, and
     method control parameters.


                               - 5 -


  -  Initialize method control parameters which may not be
     modified.

  -  Loop for each integration step.

       *  Initialize the first few series terms.

       *  Generate the entire series.

       *  Call subroutine RDCV to determine the optimal stepsize
          from (a)the location and order of the primary
          singularities, (b)the series length, (c)the error
          tolerance, and (d)adjust the stepsize.

       *  Call subroutine RSET to perform analytic continuation,
          and to print the solution.

    In ATSPGM, the stepsize used to expand the series is related to
the radius of convergence at each integration step.  After a series
is generated, the location and order of the primary singularity are
calculated.  Then, the stepsize is adjusted to control the error in
the following manner.  The terms of the series for a function g(x)
expanded at Xo with a stepsize of H := X-Xo are stored as reduced
derivatives, G(k+1) := G(k) (Xo) H**k/k!.  The stepsize H can be
varied to control the error by multiplying G(k+1) by (HNEW/H)**k.

    An exception is made in the solution of stiff problems.  The
step-size is determined in stiff problems by the length of a
polynomial that can adequately represent the function.

    A method which uses an infinite Taylor series is A-stable;
however, in practice the series must be truncated to N terms.
Then, the characteristic polynomial is p(x,y) = x - Sum[y(k)/k!].
For example, the real-valued stability intervals are (-8.85,0),
(-12.58,0), and (-16.29,0) for N = 20, 30, and 40, respectively.
Taylor series methods are best suited to solve problems with high
accuracy.  However, since very high order derivatives are used in
these methods, the solution of stiff problems can be easily solved
using the approximation of a polynomial with an exponential.




1.5 Purpose of the User Manual


    This ATOMCC User Manual is designed to support easy, and
efficient use of the ATOMCC system.  Chapter 2 may be used as a
tutorial; the rest of this User Manual is written as a reference
manual, not as a tutorial, so information is repeated as
appropriate when it applies to different issues.

    Chapter 1 presents an overview of the ATOMCC system to help you
understand how its components fit together.  This information is
helpful to using the system as described in the rest of the
manual.  A more detailed discussion can be found in references (3)
and (5).


                               - 6 -


    Chapter 2 is written as a tutorial for new users of the ATOMCC
system.  Its purpose is to show you how to use the ATOMCC system to
solve initial value problems in ODE's.  It assumes that you are
familiar with FORTRAN programming, and with the concept of comput-
ing a solution to a system of ODE's.  It gives examples showing how
to solve some specific differential equations.

    Chapter 3 is written for users who already have some experience
using the ATOMCC system.  This chapter is the heart of this User
Manual.  It attempts to show you how to use each of the features
available from the ATOMCC translator and from the ATSPGM object
program.  It is organized for reference, not for sequential
reading.




1.6 Acknowledgements


    The author would like to express his gratitude to Roy Morris
for the initial design and coding of the translator program, to
students John Fauss, David Lowery, and Manuel Prieto for their work
on series analysis, to Ray Moore, Mike Tabor, John Weiss, Mike
Ziegler, Phil Bender, and others for many helpful suggestions, and
to Jon Wright for the many suggestions which arose from his
extensive use of the ATOMCC system.  He is particularly indebted to
Professor George Corliss, who assisted him in every aspect of this
program.




1.7 References


    Here we include some references which are relevant to the
ATOMCC system.  A more complete list of references appears in
refernce(5).

  1.  D. Barton, I. M. Willers, and R. V. M. Zahar, The automatic
      solution of ordinary differential equations by the method of
      Taylor series, Comput. J., v. 14, 1971, pp. 243-248.

  2.  Y. F. Chang, Automatic solution of differential equations, in
      Constructive and Computational Methods for Differential and
      Integral Equations, edited by D. L. Colton and R. P. Gilbert,
      Lecture Notes in Math., vol. 430, Springer-Verlag, New York,
      1974, pp. 61-94.

  3.  Y. F. Chang and G. F. Corliss, Compiler for the solution of
      ordinary differential equations using Taylor series,
      Marquette University technical report, 1981.

  4.  Y. F. Chang and G. F. Corliss, Ratio-like and recurrence
      relation tests for convergence of series, J. Inst. Math.
      Appl., v. 25, 1980, pp. 349-359.


                               - 7 -


  5.  Y. F. Chang and G. F. Corliss, Solving ordinary differential
      equations using Taylor series, ACM Trans. Math. Soft, v. 8,
      1982, pp. 114-144.

  6.  Y. F. Chang, J. Fauss, M. Prieto and G. F. Corliss,
      Convergence analysis of compound Taylor series, Proceedings
      of the Seventh Conference on Numerical Mathematics and
      Computing, University of Manitoba, 1978, pp. 129-152.

  7.  Y. F. Chang, M. Tabor, J. Weiss, and G. Corliss, On the
      analytic structure of the Henon Heiles system, Phys. Lett.
      85A (1981), pp. 211-213.

  8.  G. F. Corliss, Integrating ODE's in the complex plane - Pole
      vaulting, Math. Comp., v. 35, 1980, pp. 1181-1189.

  9.  G. F. Corliss and D. Lowery, Choosing a stepsize for Taylor
      series methods for solving ODE's, J. Comput. Appl. Math., v.
      3, 1977, pp. 251-256.

  10. R. E. Moore, Interval Analysis, Prentice-Hall, Englewood
      Cliffs, NJ, 1966.

  11. L. B. Rall, Automatic Differentiation: Techniques and
      Applications, Lecture Notes in Computer Science #120,
      Springer-Verlag, Berlin, 1981.


                               - 8 -








                             Chapter 2

                           For New Users



    This Chapter is written for new users of the ATOMCC system.
Its purpose is to show how to use this system to solve initial
value problems in ordinary differential equations.  Here, we give
some specific examples of how to solve a system of differential
equations.  More detailed explanations are found in Chapter 3.




2.1 Task of the ATOMCC Translator


    The ATOMCC system is a tool to help you solve differential
equations.  It consists of two major components:- a translator
program (ATOMCC.EXE), and a subroutine object library (RDCV, DRDCV,
CRDCV, CDRDCV).  To understand the operation of these two
components, you must first understand the six steps involved in
using the system.  We discuss the purpose of each step briefly, to
acquaint you with the terms used in the detailed discussion in
Section 2.2.

    At Step 1 (edit ODEINP), the system of differential equations
are stated in the form which ATOMCC expects.  The input file ODEINP
contains four separate blocks.  (The name ODEINP is fixed within
ATOMCC; so you must use this name.)  The first block contains the
differential and algebraic equations.  ATOMCC compiler processes
the data in this block to produce a FORTRAN object program ATSPGM
which is then compiled and executed to solve the problem.  (The
name ATSPGM is also fixed within ATOMCC.)  The second, third, and
fourth blocks are copied unchanged from ODEINP to predetermined
locations in ATSPGM.

    At Step 2 (run ATOMCC), ATOMCC (a)reads the first block from
ODEINP, (b)analyzes the differential equations, and (c)copies the
second, third, and fourth blocks from ODEINP directly into ATSPGM
at locations shown by examples below.

    At Step 3 (compile ATSPGM), the ATSPGM program is compiled
using Microsoft-FORTRAN ver 3.30 compiler.

    At step 4 (link ATSPGM & RDCV), ATSPGM.OBJ is linked with the
ATOMCC subroutine library RDCV.OBJ and FORTRAN.LIB to produce an
executable module, ATSPGM.EXE.

    The recommended manner to supply the initial conditions, the
interval of integration, and control parameters is to read them


                               - 9 -


from a data file which you prepare at Step 5. The format of this
data file is completely under your control, as shown by examples
below.  Step 5 may be done at any time before Step 6, and it may be
omitted completely.

    At Step 6 (run ATSPGM), the differential equations are solved.
Each component of the equations is expanded in a Taylor series, and
the solution point is moved forward by analytic continuation.
ATSPGM reads the data file prepared at Step 5 and writes the
solution results.  The exact content, format, and location of the
solution results depend on the data in ODEINP prepared at Step 1.
Examples given below and in Chapter 3 show how this is done.

          +------------------------+
          l         Step 1         l
          l       edit ODEINP      l
          +------------------------+
                       l
          +------------------------+
          l         Step 2         l
          l     run ATOMCC.EXE     l
          +------------------------+
                       l
          +------------------------+
          l         Step 3         l       +---------------+
          l     compile ATSPGM     l       l   (RDCV.OBJ)  l
          +------------------------+       +---------------+
                       l                           l
                       l <--------------------------
                       l
          +------------------------+       +--------------+
          l         Step 4         l       l    Step 5    l
          l   link ATSPGM & RDCV   l       l prepare data l
          +------------------------+       +--------------+
                       l                           l
                       l <--------------------------
                       l
          +------------------------+
          l         Step 6         l
          l     run ATSPGM.EXE     l
          +------------------------+

                 Steps for using the ATOMCC system




2.2 Using the ATOMCC system


    In this section, we take you step-by-step through an example
using the ATOMCC system.


                              - 10 -




2.2.1 Step 1 - edit ODEINP

    The input file ODEINP specifies for ATOMCC

  1.  the system of differential equations to be solved,

  2.  how the initial conditions and the interval of integration
      are communicated to ATSPGM, and

  3.  the commands to control the operation of ATOMCC or to control
      the execution of ATSPGM.

    The statements in ODEINP follows FORTRAN conventions.  A "C" in
column 1 denotes a COMMENT, columns 1-5 are used for label numbers,
column 6 is used for continuation, columns 7-72 contain statements,
and columns 73-80 are ignored.  The statements in ODEINP may be in
either upper-case or lower-case letters.  In our discussions in
this Manual, we use upper-case letters for emphasis.  (One word of
caution:- ATOMCC does not recognize tabs.)

                 Example 2-1.  Simple ODEINP file.


C FIRST PAINLEVE TRANSCENDENT
      DIFF(Y,T,2) = 6.0*Y*Y + T    $
      $
C READ INTEGRATION INTERVAL AND INITIAL CONDITIONS.
      READ(5,1010) START,END,Y(1),Y(2)
 1010 FORMAT(4F10.3)
      WRITE(*,1020) START,END,Y(1),Y(2)
 1020 FORMAT(' SOLVE THE FIRST PAINLEVE TRANSCENDENT' /
     + ' INTERVAL:          ',2F10.3 /
     + ' INITIAL CONDITIONS:',2F10.3 /)       $
      $

    ODEINP must contain four blocks.  Each block must terminate
with the block terminator "$" in columns 7-72. Blocks 2 and 4 are
empty in Example 2-1 above.

    The first block contains the system of differential equations.
These equations are processed by ATOMCC to determine the recurrence
relations that are written into ATSPGM to generate the Taylor
series for each component of the solution.  To enter the
differential equations, DIFF(Y,X,N) is used to denote the N-th
derivative of Y with respect to X. The value of N may range from 1
to 6, inclusively.  The DIFF(,,) function is used to specify the
system of ODE's with FORTRAN-like statements using standard FORTRAN
operators and functions.

    Rarely, ATOMCC may fail to produce the correct ATSPGM for your
problem.  In such a case, write your equations differently using
many auxiliary variables.  This will allow you to solve your
problem; then, send a copy of the ODEINP that caused the problem to
Y. F. Chang, Claremont McKenna College, Claremont, CA, 91711.


                              - 11 -


    The first block can also be used to control the operation of
ATOMCC.  The most commonly used option is for ATOMCC to write
ATSPGM in double-precision with a "COPTION DOUBLE" card at the
beginning of block 1.

               Example 2-2.  Double precision ATSPGM.

COPTION DOUBLE
C
C FIRST PAINLEVE TRANSCENDENT
      DIFF(Y,T,2) = 6.0*Y*Y + T    $
      $
C READ INTEGRATION INTERVAL AND INITIAL CONDITIONS.
      READ(5,1010) START,END,Y(1),Y(2)
 1010 FORMAT(4F10.3)
      WRITE(*,1020) START,END,Y(1),Y(2)
 1020 FORMAT(' SOLVE THE FIRST PAINLEVE TRANSCENDENT' /
     + ' INTERVAL:          ',2F10.3 /
     + ' INITIAL CONDITIONS:',2F10.3 /)       $
      $

    The second block is usually empty.  It is used to insert
non-executable FORTRAN statements at the beginning of ATSPGM, such
as a SUBROUTINE card, a DIMENSION card, a COMMON card, etc.

    The second, third, and fourth blocks are not processed
syntactically by ATOMCC; they are copied directly from ODEINP into
ATSPGM.  Example 2-3 is the ATSPGM program written by ATOMCC for
the ODEINP file shown in Example 2-1. Notice where block 3 is
copied into an early part of ATSPGM.

               Example 2-3.  ATSPGM for Example 2-1.


c*+*+*+*+*+*
c     This program was produced by the  ATOMCC  translator version 7.10
c                                       Copyright (C) 1985, Y. F. Chang
c*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+
c Portions (c) Copyright, Microsoft Corp., 1981.  All rights reserved.
c This program was written for the following inputs
c
C FIRST PAINLEVE TRANSCENDENT
C     DIFF(Y,T,2) = 6.0*Y*Y + T
c--------
c no instructions in second input block
c--------
      COMMON /IPASS/ LENSER,LENVAR,MPRINT,MSTIFF,LRUN,
     + KTRDCV,KNTSTP,KTSTIF,KXPNUM,KDIGS,KENDFG,NTERMS,NOPT
     A /RPASS/ RADIUS,ERRLIM,ADJSTF,RCREAL,RCIMAG,RDCERR
     B /CPASS/ START,END,ORDER
     C /DPASS/ H,HNEW,XPRINT,DLTXPT
      DIMENSION TMPS(  36,  1)
      CHARACTER*6 NAMES
      EQUIVALENCE (TMPS(1,1),Y(1))
      DIMENSION NAMES(1), Y(36), T(2), TMPAAB(30), TMPAAA(30)
      DATA NAMES(1)/'Y.....'/
      Y(33) = 1.1
   10 FORMAT(72H   ATOMCC  Ver. 7.10, Copyright (C) 1985, Y. F. Chang; S


                              - 12 -


     Aolution results./9H   ******)
   11 FORMAT(/5X,11HStep number,I6,13H at the point,1P1E12.4/1X,
     A 9Hvalues of )
   12 FORMAT(1X, A6,(1X,1P4E13. 5))
   13 FORMAT(5X,21HStepsize adjusted to ,1PE13.5)
   14 FORMAT(/5X,35HThe solution stopped normally after, I4,24H steps as
     a set by nsteps. )
   16 FORMAT(/5X,63HThe adjustment for stepsize seems to be in a loop. P
     Alease try a /5X,22Hshorter series length. )
      WRITE(*,10)
c--------
c Initialize variables to default values.
c--------
      NSTEPS = 40
      H = 1.E0
      ERRLIM = 1.E- 6
      LENSER = 30
      MPRINT = 4
      NTERMS = 2
      KTRDCV =  1
      ADJSTF = 1.E-2
      MSTIFF =  0
      DLTXPT = 0.E0
c--------
c start of third  input block
c--------
C READ INTEGRATION INTERVAL AND INITIAL CONDITIONS.
      READ(5,1010) START,END,Y(1),Y(2)
 1010 FORMAT(4F10.3)
      WRITE(*,1020) START,END,Y(1),Y(2)
 1020 FORMAT(' SOLVE THE FIRST PAINLEVE TRANSCENDENT' /
     + ' INTERVAL:          ',2F10.3 /
     + ' INITIAL CONDITIONS:',2F10.3 /)
c--------
c end of third  input block
c--------
c More initializations
c--------
      DLTXPT = SIGN(DLTXPT,(END-START))
      H = SIGN(H,(END-START))
      KDIGS =  6
      XPRINT = START + DLTXPT
      KXPNUM = 35
      LENVAR = 36
      LRUN = 1
      KTSTIF = 0
      NUMEQS =   1
      IF(LENSER.GT.(LENVAR- 6)) LENSER = LENVAR -  6
      IF(MPRINT.LT.2) GO TO 17
      WRITE(*,11) KTSTIF,START
         K = Y(33)
         WRITE(*,12) NAMES(K),Y(1), Y(2)
c--------
c Loop for integration steps.  Inside the loop, print the desired output
c--------
   17 DO 27 KINTS=1,NSTEPS
         KOUNT = 0
         KNTSTP = KINTS


                              - 13 -


   19    CONTINUE
         T(1) = START
         T(2) = H
         Y(2) = Y(2)*(H)
c--------
c Preliminary series calculations
c--------
         TMPAAA(1) = 6.E0*Y(1)
         TMPAAB(1) = TMPAAA(1)*Y(1)
         Y(3) = (TMPAAB(1) + T(1))*(H*H/2.E0)
         TMPAAA(2) = 6.E0*Y(2)
         TMPAAB(2) = TMPAAA(1)*Y(2) + TMPAAA(2)*Y(1)
         Y(4) = (TMPAAB(2) + T(2))*(H*H/6.E0)
c--------
c Loop for series calculations
c--------
         DO 23 K= 5,LENSER
            KA = K - 1
            KB = K - 2
            TMPAAA(KB) = 6.E0*Y(KB)
            TMPAAB(KB) = 0.E0
            KZ = 1 + KB
            DO 30 N=1, KB
               L = KZ - N
               TMPAAB(KB) = TMPAAB(KB) + TMPAAA(N)*Y(L)
   30            CONTINUE
            Y(K) = (TMPAAB(KB))*(H*H/(KB*KA))
c--------
c Test and adjust H to avoid over/under flow.
c--------
            IF(MSTIFF.GE.20 .AND. KTSTIF.GT.0) GO TO 23
            TMP = ABS(Y(K))
            IF(TMP.LE.1.E-35) GO TO 23
            IF(TMP.LT.1.E20 .AND. TMP.GT.1.E-20) GO TO 23
            IF(KTSTIF.NE.0 .AND. TMP.LT.1.0) GO TO 23
            KOUNT = KOUNT + 1
            IF(KOUNT.LT.9) GO TO 22
            WRITE(*,16)
            GO TO 28
   22       CONTINUE
            Y(2) = Y(2)/(H)
            H = H * TMP**(0.3/(1-K))
            IF(MPRINT.GE.4) WRITE(*,13) H
            GO TO 19
   23    LRUN = 1
c--------
c Calculate radius of convergence and take optimum step.
c--------
         CALL RDCV(TMPS,LENVAR,NUMEQS,NAMES)
   24    CALL RSET(TMPS,LENVAR,NUMEQS,NAMES)
c--------
c no instructions in fourth input block
c--------
   25    GO TO (26,28,24), KENDFG
   26    H = SIGN(RADIUS,H)
         START = START + HNEW
         IF(MPRINT.LT.4) GO TO 27
         WRITE(*,11) KNTSTP, START


                              - 14 -


         K = Y(33)
         WRITE(*,12) NAMES(K),Y(1), Y(2)
   27 CONTINUE
      WRITE(*,14) NSTEPS
   28 CONTINUE
   29 STOP
      END

    The third block is usually used to specify the interval of
integration and the initial conditions by reading them from a data
file prepared at Step 5. This is the file DATA opened in block 3.
The interval of integration is from START to END.  END is allowed
to be less than START for integration in a negative direction.  The
initial values (at START) of a dependent variable named y and its
derivatives are assigned to the array Y as follows:-

      Y(1) denotes y at START,
      Y(2) denotes y' at START,
      Y(3) denotes y'' at START, etc.

Thus in Example 2-1, two initial conditions Y(1) for y(0) and Y(2)
for y'(0) are entered for the second order differential equation.

    Any valid FORTRAN statement may be included in block 3 to be
copied into ATSPGM, as shown in Example 2-1 by the WRITE statement
to echo the input.  The third block may also be used to change the
default values of method-controlling variables.  You can see in
Example 2-3 that many variables are initialized immediately before
block 3. The meaning of these variables is described in Chapter 3.

    The fourth block is usually empty.  It may be used to insert
statements into ATSPGM at the end of each integration step.

    This concludes the discussion of how to prepare the input
file.  More information about the use of specific features can be
found in Chapter 3.


2.2.2 Step 2 - Run ATOMCC

    The appropriate command to execute the ATOMCC compiler is
simply [ATOMCC].  The ATOMCC translator uses two files:- ODEINP for
the equation statements and initial conditions, and ATSPGM for the
output object program.  (The names ODEINP and ATSPGM are fixed
within ATOMCC.)  The messages produced by ATOMCC are placed on your
terminal.  To have the messages written onto a disc file (say MSG),
use the command [ATOMCC > MSG].


                              - 15 -



         Example 2-4.  Translator messages for Example 2-1.


              ATOMCC Ver. 7.10, Copyright (C) 1985, Y. F. Chang.
              --------------------------------------------------
 Portions (c) Copyright, Microsoft Corp., 1981.  All rights reserved.


  FIRST PAINLEVE TRANSCENDENT
      DIFF(Y,T,2) = 6.0*Y*Y + T    $
  equation   1 is in position   1
      $
  READ INTEGRATION INTERVAL AND INITIAL CONDITIONS.
      READ(5,1010) START,END,Y(1),Y(2)
 1010 FORMAT(4F10.3)
      WRITE(*,1020) START,END,Y(1),Y(2)
 1020 FORMAT(' SOLVE THE FIRST PAINLEVE TRANSCENDENT' /
     + ' INTERVAL:          ',2F10.3 /
     + ' INITIAL CONDITIONS:',2F10.3 /)       $
      $
          ATOMCC completed
 Stop - Program terminated.


2.2.3 Step 3 and 4 - Compile and link ATSPGM

    As you get comfortable with the ATOMCC system, you will rarely
inspect ATSPGM, unless either an error occurs, or you choose to
edit ATSPGM by hand.

    ATSPGM, written by ATOMCC, is just like any other FORTRAN
program; you may edit it to suit your needs.  Whether edited or
not, ATSPGM is ready to be compiled and linked with the necessary
subroutines from the ATOMCC library (RDCV).

    The appropriate commands to compile and link ATSPGM are:-

  -  B:FOR1 ATSPGM.  ;

  -  B:PAS2

  -  B:LINK ATSPGM+RDCV,,NUL,B:


2.2.4 Step 5 - Prepare the data

    At Step 1, when you prepared ODEINP for ATOMCC, you may have
included some READ statements in block 3 to communicate the
interval of integration and the initial conditions to ATSPGM.
Before you run ATSPGM, the data file to be read by those statements
must be prepared with the appropriate file name given in your OPEN
statement.


                              - 16 -



              Example 2-5.  Data file for Example 2-1.

     0.000     1.100     1.000     0.000

    If you know that you will be solving a simple problem only
once, Step 5 can be eliminated by stating the values of START, END,
and the initial conditions with FORTRAN assignment statements in
block 3 as shown below.

          Example 2-6.  Assignment statements in block 3.


C First Painleve transcendent
      DIFF(Y,T,2) = 6.0*Y*Y + T    $
      $
C Assign integration interval and initial conditions.
      START = 0.0
      END = 1.1
      Y(1) = 1.0
      Y(2) = 0.0
      WRITE(*,1020) START,END,Y(1),Y(2)
 1020 FORMAT(' Solve the first Painleve transcendent' /
     +  ' Interval:          ',2F10.3 /
     +  ' Initial conditions:',2F10.3 /)       $
      $


2.2.5 Step 6 - Run ATSPGM

    At step 6, you are ready to run ATSPGM; the command is simply
[ATSPGM].  ATSPGM writes its output to your terminal, for solution
output on a disc file (say PRTOUT) use [ATSPGM > PRTOUT].

           Example 2-7.  Solution Output for Example 2-1.


  ATOMCC  Ver. 7.10, Copyright (C) 1985, Y. F. Chang; Solution results.
  ******
SOLVE THE FIRST PAINLEVE TRANSCENDENT
INTERVAL:                .000     1.100
INITIAL CONDITIONS:     1.000      .000

    Step number     0 at the point   .0000E+00
values of
Y.....   1.00000E+00   .00000E+00

    Step number     1 at the point  7.1000E-01
values of
Y.....   4.04877E+00  1.62701E+01


                              - 17 -




    Step number     2 at the point  1.0100E+00
values of
Y.....   2.58324E+01  2.62656E+02

    Step number     3 at the point  1.1000E+00
values of
Y.....   8.77693E+01  1.64459E+03
Stop - Program terminated.




2.3 Output at equally spaced points


    You should be able to use ATOMCC to solve routine problems.
The points at which ATSPGM computes the solutions are determined by
the actual integration steps taken, which are not uniform in size.
This Section shows you how to force ATSPGM to print the solutions
at equally spaced points.

    The output from ATSPGM is controlled by two variables, MPRINT
(amount of print), and DLTXPT (print interval).  To produce output
at equally spaced points, assign MPRINT=2 to turn off the print at
the actual integration steps, and assign DLTXPT=DELTA, where DELTA
is your desired print interval.

            Example 2-8.  Equally spaced output points.


c First Painleve transcendent
      DIFF(Y,T,2) = 6.0*Y*Y + T    $
      $
c Assign integration interval and initial conditions.
      MPRINT = 2
      DLTXPT = 0.2
      START = 0.0
      END = 1.1
      Y(1) = 1.0
      Y(2) = 0.0
      WRITE(*,1020) START,END,Y(1),Y(2)
 1020 FORMAT(' Solve the first Painleve transcendent' /
     +  ' Interval:          ',2F10.3 /
     +  ' Initial conditions:',2F10.3 /)       $
      $


2.3.1 ZEROT - Stopping at roots of variables

    It is often of interest to locate points at which a component
of the solution has a root or assumes some specified value.  The
subroutine ZEROT automatically solve such problems.  DZEROT is the
double-precision version.


                              - 18 -


    The form of the CALL is

             CALL ZEROT(NUMBER,Y,ROOT,KEY,TMPS,LENVAR,NUMEQS)

    where

    NUMBER is the index of the Y series term whose root is sought,
    Y      is the variable whose root is desired,
    ROOT   is the value Y is to assume (= 0 for a root),
    KEY    is 1 if Y is a dependent variable, or
              0 if Y is not a dependent variable.

    The arguements TMPS, LENVAR, and NUMEQS must be exactly as
written above.

               Example 3-16.  Rootfinding with ZEROT.


      DIFF(Y,T,2) = 6*Y*Y + T  $
      $
      START = 0.0
      END = 1.15
      ROOT = 20.0
      Y(1) = 1.0
      Y(2) = 0.0
      $
      CALL ZEROT(1,Y,ROOT,1)
      $

    When the variable whose root is being sought is not a dependent
variable, KEY is set to 0.

        CALL ZEROT(2,VARY,0.0,0)
        IF(LRUN.NE.0) GO TO 25
        TEMP = START + HNEW
        WRITE(7,1010) KINTS,TEMP,VARY(1),VARY(2)
   1010 FORMAT(I5,3F10.4)
        GO TO 25  $

    In these examples, it is not necessary for one to print the
information as shown in the second case.  The ATSPGM program does
stop and restart the solution automatically at the exact root and
the output is controlled by MPRINT.

    The index NUMBER can be set at any positive (non-zero) integer
value; however, obviously when NUMBER is very large the accuracy of
the root will suffer.


2.3.2 MSTIFF=20,21,22 - Stiff problems.

    This version of ATOMCC contains a double-precision algorithm to
solve stiff problems.  To use it, one can either set MSTIFF=20, or
21, or 22. Other parameters that should be controlled are H,
ADJSTF, and NSTEPS.  It is also desirable to set MPRINT to 7, at
least initially, for observing the progress of the solution.  If it
should be evident that the problem is not really stiff, then it is
most advisable to solve it as a normal problem.


                              - 19 -


    MSTIFF=20 is the more conservative of the three algorithms.  In
this case, LENSER is set to be 15. The default value for ADJSTF,
the error-controlling parameter, is a rather large 1.E-2. The user
should run the stiff solution at least one more time with a
somewhat smaller ADJSTF, say 1.E-3, to check on its validity.

    When MSTIFF=21, LENSER is set to only 10. So, this option
should be used only if the user is absolutely certain that the
problem under study is very stiff.  The solution of stiff problems
under this option is considerably faster than that for MSTIFF=20.

    MSTIFF=22 is identical to MSTIFF=20 except for the fact that
there is no attempt to identify steady-state solutions.

    As mentioned above, the stiff algorithm is written in double-
precision.  It is simply not cost effective to solve such problems
using single-precision.  There is one other restriction on stiff
problems.  All such problems must be stated as first-degree ODE's.


                              - 20 -








                             Chapter 3

                          How to Use ----



    This Chapter is written for users who already have solved
several problems using the ATOMCC system.  It assumes familiarity
with FORTRAN programming, with Chapter 2, and with the numerical
solution of ODE's.  This Chapter is the heart of this User Manual.
It attempts to show how to use each of the features available from
ATOMCC and from ATSPGM.  It is organized for reference, not for
sequential reading.  Consequently, some information found in other
parts of this Manual are repeated here.




3.1 Solving your problem


    The tasks which must be accomplished in order to run the ATOMCC
system on your computer were discussed in Section 2.2. They are:-

  Edit     ODEINP      (containing ODE's)

  Run      ATOMCC.EXE  (execute ATOMCC translator)
           (This creates ATSPGM, the object FORTRAN program,
            which is treated like any FORTRAN program.
            ATSPGM may be edited.)

  Compile  ATSPGM.
           (This creates ATSPGM.OBJ, the object module.)

  Link     ATSPGM.OBJ,  with library options
                        (RDCV.OBJ    for single precision
                         DRDCV.OBJ   for double precision
                         CRDCV.OBJ   for complex, and
                         CDRDCV.OBJ  for complex double)
           (This creates ATSPGM.EXE, the execution module.)

  Edit     DATA input-file if any is used.

  Run      ATSPGM.EXE


3.1.1 Translator file, ODEINP

    The ODEINP file contains the ODE's to be solved and information
specifying how the initial conditions and the integration interval
are determined in ATSPGM.  It may contain commands to control
(a)the operation of the ATOMCC translator, (b)the execution of the


                              - 21 -


solution, and (c)the desired format of the output.  The names
ODEINP and ATSPGM are fixed within ATOMCC; you must have these
files on your disc when you execute ATOMCC.

    The data in ODEINP follows FORTRAN conventions.  Columns 1-5
are used for line numbers, column 6 is used for continuation
characters, columns 7-72 contain statements, and columns 73-80 are
ignored.  As in FORTRAN, all blanks are ignored.  A 'C' in column 1
denotes a comment which is copied directly into the ATSPGM file.
The statements in ODEINP can be either upper-case or lower-case
letters.  We use upper-case in this manual for emphasis.  (A word
of caution, the tab character is not recognized by ATOMCC.)

    The ODEINP file contains four blocks.  Each block ends with the
block terminator symbol '$' in columns 7-72. (A comment card must
not contain a block terminator '$'.)  Sections 3.2-3.4 discuss each
of the blocks in detail.  Here is an example of an input file which
illustrates several of the features which will be discussed in
Sections 3.2-3.5.

                     Example 3-1.  ODEINP file.


C Block 1
C
C  System with parameter.
C
      DIFF(X,T,2) = - ALPHA*X*R
      DIFF(Y,T,2) = - ALPHA*Y*R
      R = (X*X + Y*Y)**(-1.5)
      ALPHA = 0.65            $
c
c Block 2
c
      CHARACTER*80 LINE       $
c
c Block 3
c
c  Read:- heading line, print code, maximum number of integration
c         steps.
c  Echo the above.
c  Read:- heading line, integration interval, print interval,
c         parameter in equations, initial conditions.
c  Echo the above.
c
      OPEN(5,FILE='DATA')
      OPEN(7,FILE='PLOTS',STATUS='NEW')
      READ(5,1010) LINE,MPRINT,NSTEPS
 1010 FORMAT(A80/2I10)
      WRITE(*,1010) LINE,MPRINT,NSTEPS
      READ(5,1020) LINE,START,END,DLTXPT,ALPHA,X(1),X(2),Y(1),Y(2)
 1020 FORMAT(A80/8F10.3)
      WRITE(*,1020) LINE,START,END,DLTXPT,ALPHA,X(1),X(2),Y(1),Y(2)
c  Assignment statements for the error control parameter
      ERRLIM = 1.0E-04        $


                              - 22 -



c
c Block 4
c
c  Produce file of data for plotting.
c
      IF(KENDFG.EQ.3) WRITE(7,1030) KINTS,XPRINT,X(1),X(2),Y(1),Y(2)
 1030 FORMAT(I5,1P5E14.5)     $

    If you have a simple problem which will be solved only once,
the contents of the 'DATA' file may be entered directly as data
into block 3. However, the compilation and linking of the ATSPGM
file takes an appreciable amount of time and therefore should be
avoided on problems that is solved more than once.

    There are two ways to solve a given problem containing
functions.  Either they can be placed as FORTRAN statements in
ODEINP to be copied directly into ATSPGM, or they can be inserted
by editing ATSPGM.  The choice depends on your personal style; the
authors prefer to work with ODEINP, because it is short.
Therefore, it is easy to find the correct place to make changes.
The cost of re-running the ATOMCC translator is well worth the
convenience, because ATOMCC is very fast.


3.1.2 Translator file, the terminal

    The translator messages contains the information which the
ATOMCC expects you to inspect.  It includes an echo of the input
file and any error messages.  The Appendix contains a list of the
error messages which may be produced.  The messages will appear on
your terminal.

         Example 3-2.  Translator messages for Example 3-1.


             ATOMCC Ver. 7.10, Copyright (C) 1985, Y. F. Chang.
             --------------------------------------------------
Portions (c) Copyright, Microsoft Corp., 1981.  All rights reserved.


c Block 1

c  System with parameter.

     DIFF(X,T,2) = - ALPHA*X*R
     DIFF(Y,T,2) = - ALPHA*Y*R
     R = (X*X + Y*Y)**(-1.5)
     ALPHA = 0.65            $
 equation   3 is in position   1
 equation   4 is in position   2
 equation   1 is in position   3
 equation   2 is in position   4

c Block 2

     CHARACTER*80 LINE       $


                              - 23 -


c Block 3

c  Read:- heading line, print code, maximum number of integration
c         steps.
c  Echo the above.
c  Read:- heading line, integration interval, print interval,
c         parameter in equations, initial conditions.
c  Echo the above.

     OPEN(5,FILE='DATA')
     OPEN(7,FILE='PLOTS',STATUS='NEW')
     READ(5,1010) LINE,MPRINT,NSTEPS
1010 FORMAT(A80/2I10)
     WRITE(*,1010) LINE,MPRINT,NSTEPS
     READ(5,1020) LINE,START,END,DLTXPT,ALPHA,X(1),X(2),Y(1),Y(2)
1020 FORMAT(A80/8F10.3)
     WRITE(*,1020) LINE,START,END,DLTXPT,ALPHA,X(1),X(2),Y(1),Y(2)
c  Assignment statements for the error control parameter
     ERRLIM = 1.0E-04        $

c Block 4

c  Produce file of data for plotting.

     IF(KENDFG.EQ.3) WRITE(7,1030)KINTS,XPRINT,X(31),X(32),Y(31),Y(32)
1030 FORMAT(I5,1P5E14.5)     $
         ATOMCC completed
Stop - Program terminated.


3.1.3 Translator file, ATSPGM

    The ATSPGM file contains the FORTRAN object program written by
ATOMCC to solve the system of differential equations using long
Taylor series.  ATOMCC uses the variable names given in ODEINP, so
that ATSPGM appears to have been custom written for the specific
problem.  Usually you do not need to inspect ATSPGM, but sometimes
you may find it necessary to edit it like you would edit any other
FORTRAN program to achieve some particular result.  Section 3.6
contains an example for which editing of ATSPGM is necessary.

               Example 3-3.  ATSPGM for Example 3-1.


c*+*+*+*+*+*
c     This program was produced by the  ATOMCC  translator version 7.10
c                                       Copyright (C) 1985, Y. F. Chang
c*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+
c Portions (c) Copyright, Microsoft Corp., 1981.  All rights reserved.
c This program was written for the following inputs
c
C BLOCK 1
C  SYSTEM WITH PARAMETER.
C     DIFF(X,T,2) = - ALPHA*X*R
C     DIFF(Y,T,2) = - ALPHA*Y*R
C     R = (X*X + Y*Y)**(-1.5)
C     ALPHA = 0.65
c--------


                              - 24 -


c start of second input block
c--------
c
c Block 2
c
      CHARACTER*80 LINE
c--------
c end of second input block
c--------
      COMMON /IPASS/ LENSER,LENVAR,MPRINT,MSTIFF,LRUN,
     + KTRDCV,KNTSTP,KTSTIF,KXPNUM,KDIGS,KENDFG,NTERMS,NOPT
     A /RPASS/ RADIUS,ERRLIM,ADJSTF,RCREAL,RCIMAG,RDCERR
     B /CPASS/ START,END,ORDER
     C /DPASS/ H,HNEW,XPRINT,DLTXPT
      DIMENSION TMPS(  36,  2)
      CHARACTER*6 NAMES
      EQUIVALENCE (TMPS(1,1),X(1)),(TMPS(1,2),Y(1))
      DIMENSION NAMES(2), T(2), Y(36), X(36), R(30), TMPAAH(30),
     A  TMPAAG(30), TMPAAF(30), TMPAAE(30), TMPAAD(30), TMPAAC(30),
     B  TMPAAB(30)
      DATA NAMES(1)/'X.....'/
      DATA NAMES(2)/'Y.....'/
      X(33) = 1.1
      Y(33) = 2.1
   10 FORMAT(72H   ATOMCC  Ver. 7.10, Copyright (C) 1985, Y. F. Chang; S
     Aolution results./9H   ******)
   11 FORMAT(/5X,11HStep number,I6,13H at the point,1P1E12.4/1X,
     A 9Hvalues of )
   12 FORMAT(1X, A6,(1X,1P4E13. 5))
   13 FORMAT(5X,21HStepsize adjusted to ,1PE13.5)
   14 FORMAT(/5X,35HThe solution stopped normally after, I4,24H steps as
     a set by nsteps. )
   16 FORMAT(/5X,63HThe adjustment for stepsize seems to be in a loop. P
     Alease try a /5X,22Hshorter series length. )
      WRITE(*,10)
c--------
c Initialize variables to default values.
c--------
      NSTEPS = 40
      H = 1.E0
      ERRLIM = 1.E- 6
      LENSER = 30
      MPRINT = 4
      NTERMS = 2
      KTRDCV =  2
      ADJSTF = 1.E-2
      MSTIFF =  0
      DLTXPT = 0.E0
c--------
c constant expressions
c--------
      ALPHA = 6.5E-1
c--------
c start of third  input block
c--------
c
c Block 3
c


                              - 25 -


c  Read:- heading line, print code, maximum number of integration
c         steps.
c  Echo the above.
c  Read:- heading line, integration interval, print interval,
c         parameter in equations, initial conditions.
c  Echo the above.
c
      OPEN(5,FILE='DATA')
      OPEN(7,FILE='PLOTS',STATUS='NEW')
      READ(5,1010) LINE,MPRINT,NSTEPS
 1010 FORMAT(A80/2I10)
      WRITE(*,1010) LINE,MPRINT,NSTEPS
      READ(5,1020) LINE,START,END,DLTXPT,ALPHA,X(1),X(2),Y(1),Y(2)
 1020 FORMAT(A80/8F10.3)
      WRITE(*,1020) LINE,START,END,DLTXPT,ALPHA,X(1),X(2),Y(1),Y(2)
c  Assignment statements for the error control parameter
      ERRLIM = 1.0E-04
c--------
c end of third  input block
c--------
      TMPAAA = -1.5E0
c--------
c constant expressions
c--------
c More initializations
c--------
      DLTXPT = SIGN(DLTXPT,(END-START))
      H = SIGN(H,(END-START))
      KDIGS =  6
      XPRINT = START + DLTXPT
      KXPNUM = 35
      LENVAR = 36
      LRUN = 1
      KTSTIF = 0
      NUMEQS =   2
      IF(LENSER.GT.(LENVAR- 6)) LENSER = LENVAR -  6
      IF(MPRINT.LT.2) GO TO 17
      WRITE(*,11) KTSTIF,START
         K = X(33)
         WRITE(*,12) NAMES(K),X(1), X(2)
         K = Y(33)
         WRITE(*,12) NAMES(K),Y(1), Y(2)
c--------
c Loop for integration steps.  Inside the loop, print the desired output
c--------
   17 DO 27 KINTS=1,NSTEPS
         KOUNT = 0
         KNTSTP = KINTS
   19    CONTINUE
         T(1) = START
         T(2) = H
         X(2) = X(2)*(H)
         Y(2) = Y(2)*(H)
c--------
c Preliminary series calculations
c--------
         TMPAAB(1) = X(1)*X(1)
         TMPAAC(1) = Y(1)*Y(1)


                              - 26 -


         TMPAAE(1) = ALPHA*X(1)
         TMPAAG(1) = ALPHA*Y(1)
         TMPAAD(1) = TMPAAB(1) + TMPAAC(1)
         R(1) = TMPAAD(1) ** TMPAAA
         TMPAAF(1) = TMPAAE(1)*R(1)
         TMPAAH(1) = TMPAAG(1)*R(1)
         X(3) = (-TMPAAF(1))*(H*H/2.E0)
         Y(3) = (-TMPAAH(1))*(H*H/2.E0)
         TMPAAB(2) = X(1)*X(2) + X(2)*X(1)
         TMPAAC(2) = Y(1)*Y(2) + Y(2)*Y(1)
         TMPAAE(2) = ALPHA*X(2)
         TMPAAG(2) = ALPHA*Y(2)
         TMPAAD(2) = TMPAAB(2) + TMPAAC(2)
         R(2) = TMPAAA*R(1)*TMPAAD(2)/TMPAAD(1)
         TMPAAF(2) = TMPAAE(1)*R(2) + TMPAAE(2)*R(1)
         TMPAAH(2) = TMPAAG(1)*R(2) + TMPAAG(2)*R(1)
         X(4) = (-TMPAAF(2))*(H*H/6.E0)
         Y(4) = (-TMPAAH(2))*(H*H/6.E0)
c--------
c Loop for series calculations
c--------
         DO 23 K= 5,LENSER
            KA = K - 1
            KB = K - 2
            KC = K - 3
            TMPAAB(KB) = 0.E0
            TMPAAC(KB) = 0.E0
            KZ = 1 + KB
            DO 30 N=1, KB
               L = KZ - N
               TMPAAB(KB) = TMPAAB(KB) + X(N)*X(L)
               TMPAAC(KB) = TMPAAC(KB) + Y(N)*Y(L)
   30            CONTINUE
            TMPAAE(KB) = ALPHA*X(KB)
            TMPAAG(KB) = ALPHA*Y(KB)
            TMPAAD(KB) = TMPAAB(KB) + TMPAAC(KB)
            R(KB) = R(1)*TMPAAD(KC+1)*(KC)*TMPAAA
            KY = 2 + KC
            DO 31 N=2, KC
               L = KY - N
               AL = (L - 1)
   31       R(KB) = R(KB) + R(N)*TMPAAD(L)*AL
     A *TMPAAA - TMPAAD(N)*R(L)*AL
            R(KB) = R(KB) /(KC)/TMPAAD(1)
            TMPAAF(KB) = 0.E0
            TMPAAH(KB) = 0.E0
            KZ = 1 + KB
            DO 32 N=1, KB
               L = KZ - N
               TMPAAF(KB) = TMPAAF(KB) + TMPAAE(N)*R(L)
               TMPAAH(KB) = TMPAAH(KB) + TMPAAG(N)*R(L)
   32            CONTINUE
            X(K) = (-TMPAAF(KB))*(H*H/(KB*KA))
            Y(K) = (-TMPAAH(KB))*(H*H/(KB*KA))
c--------
c Test and adjust H to avoid over/under flow.
c--------
            IF(MSTIFF.GE.20 .AND. KTSTIF.GT.0) GO TO 23


                              - 27 -


            TMP = ABS(X(K))
            IF(TMP.LE.1.E-35) GO TO 23
            IF(TMP.LT.1.E20 .AND. TMP.GT.1.E-20) GO TO 23
            IF(KTSTIF.NE.0 .AND. TMP.LT.1.0) GO TO 23
            KOUNT = KOUNT + 1
            IF(KOUNT.LT.9) GO TO 22
            WRITE(*,16)
            GO TO 28
   22       CONTINUE
            X(2) = X(2)/(H)
            Y(2) = Y(2)/(H)
            H = H * TMP**(0.3/(1-K))
            IF(MPRINT.GE.4) WRITE(*,13) H
            GO TO 19
   23    LRUN = 1
c--------
c Calculate radius of convergence and take optimum step.
c--------
         CALL RDCV(TMPS,LENVAR,NUMEQS,NAMES)
   24    CALL RSET(TMPS,LENVAR,NUMEQS,NAMES)
c--------
c start of fourth input block
c--------
c
c Block 4
c
c  Produce file of data for plotting.
c
      IF(KENDFG.EQ.3) WRITE(7,1030)KINTS,XPRINT,X(31),X(32),Y(31),Y(32)
 1030 FORMAT(I5,1P5E14.5)
c--------
c end of fourth input block
c--------
   25    GO TO (26,28,24), KENDFG
   26    H = SIGN(RADIUS,H)
         START = START + HNEW
         IF(MPRINT.LT.4) GO TO 27
         WRITE(*,11) KNTSTP, START
         K = X(33)
         WRITE(*,12) NAMES(K),X(1), X(2)
         K = Y(33)
         WRITE(*,12) NAMES(K),Y(1), Y(2)
   27 CONTINUE
      WRITE(*,14) NSTEPS
   28 CONTINUE
   29 STOP
      END


3.1.4 DATA input file

    For most problems solved using the ATOMCC system, you should
communicate the initial conditions, the integration interval, the
coefficients in the differential equations, and the ATOMCC control
parameters by reading them from a DATA input file.


                              - 28 -



           Example 3-4.  DATA input file for Example 3-1.


  MPRINT   NSTEPS           for example 3-1
      4       80
  START      END   DLTXPT   ALPHA     X(1)     X(2)     Y(1)     Y(2)
  1.0     10.0      0.25     0.58    -1.0      0.0      0.0      4.3


3.1.5 Solution file

    ATSPGM writes all of its messages and answers to your
terminal.  The format of the solution depends on the ATOMCC control
parameters (see Section 3.4). If you should wish to have the
solution written onto a disc file, you can use the execution
statement (ATSPGM > PRTOUT).  Then, all messages and answers will
be written to the solution file PRTOUT.  The solution file can also
contain informa- tion written by user supplied WRITE(*,xxx)
statements.  A portion of the solution file is given below.

            Example 3-5.  Solution file for Example 3-1.


  ATOMCC  Ver. 7.10, Copyright (C) 1985, Y. F. Chang; Solution results.
  ******
 MPRINT   NSTEPS           for example 3-1
     4       80
 START      END   DLTXPT   ALPHA     X(1)     X(2)     Y(1)     Y(2)
 1.000   10.000    .250     .580   -1.000     .000     .000    4.300

    Step number     0 at the point  1.0000E+00
values of
X.....  -1.00000E+00   .00000E+00
Y.....    .00000E+00  4.30000E+00

    Step number     1 at the point  1.1430E+00
values of
X.....  -9.94536E-01  7.08452E-02
Y.....   6.13850E-01  4.27990E+00

      Step number     2 at the point  1.2500E+00
  values of
  X.....    -9.85268E-01  9.92472E-02
  Y.....     1.07052E+00  4.25646E+00

    Step number     2 at the point  1.3400E+00
values of
X.....  -9.75709E-01  1.11976E-01
Y.....   1.45285E+00  4.24032E+00

      Step number     3 at the point  1.5000E+00
  values of
  X.....    -9.56779E-01  1.23037E-01
  Y.....     2.12958E+00  4.22039E+00


                              - 29 -




    Step number     3 at the point  1.6400E+00
values of
X.....  -9.39208E-01  1.27496E-01
Y.....   2.71959E+00  4.20915E+00

      Step number     4 at the point  1.7500E+00
  values of
  X.....    -9.25064E-01  1.29523E-01
  Y.....     3.18223E+00  4.20277E+00

      Step number     4 at the point  2.0000E+00
  values of
  X.....    -8.92333E-01  1.31982E-01
  Y.....     4.23159E+00  4.19295E+00

    Step number     4 at the point  2.1300E+00
values of
X.....  -8.75128E-01  1.32677E-01
Y.....   4.77644E+00  4.18942E+00


3.1.6 User files

    It is often useful to create your own output files using WRITE
statements to produce output in a format of your choice.  Example
3-6 shows a portion of such a file for plotting, produced by
Example 3-1.

               Example 3-6.  User file for plotting.


2   1.25000E+00  -9.85268E-01   9.92472E-02   1.07052E+00   4.25646E+00
3   1.50000E+00  -9.56779E-01   1.23037E-01   2.12958E+00   4.22039E+00
4   1.75000E+00  -9.25064E-01   1.29523E-01   3.18223E+00   4.20277E+00
4   2.00000E+00  -8.92333E-01   1.31982E-01   4.23159E+00   4.19295E+00
5   2.25000E+00  -8.59178E-01   1.33133E-01   5.27900E+00   4.18678E+00
5   2.50000E+00  -8.25810E-01   1.33750E-01   6.32514E+00   4.18258E+00
6   2.75000E+00  -7.92324E-01   1.34112E-01   7.37039E+00   4.17953E+00
6   3.00000E+00  -7.58765E-01   1.34340E-01   8.41497E+00   4.17723E+00
6   3.25000E+00  -7.25160E-01   1.34490E-01   9.45904E+00   4.17542E+00
6   3.50000E+00  -6.91524E-01   1.34594E-01   1.05027E+01   4.17398E+00
6   3.75000E+00  -6.57866E-01   1.34667E-01   1.15461E+01   4.17279E+00
7   4.00000E+00  -6.24192E-01   1.34720E-01   1.25891E+01   4.17179E+00




3.2 Using block 1


    The first block contains the system of differential equations.
These equations are processed by ATOMCC to determine the recurrence
relations which should be written into ATSPGM to generate the
Taylor series for each component of the solution.


                              - 30 -


    The first block may also be used for ATOMCC options to control
the operation of the translator.  This is done by placing a COPTION
card at the beginning of ODEINP.  Multiple translator options may
be specified on the same line, i.e. COPTION DOUBLE, LENVAR=40.
Multiple COPTION cards are also allowed, but they must precede all
other cards.


3.2.1 Format for the system of equations

    For the differential equations, DIFF(Y,X,N) is used to denote
the n-th derivative of the dependent variable y with respect to the
independent variable x.  The value of the variable N may range from
1 to 6, inclusively.  The DIFF(X,Y,N) function is used to specify
ODE's just as with other standard FORTRAN operators and functions.
All functions supported by ATOMCC are listed in Section 4.2. The
statements in ODEINP can be either upper- or lower-case letters.
We use upper-case in this manual for emphasis.

    The input to ATOMCC follows FORTRAN conventions.  Comment cards
contain a 'C' in column 1. The entire comment card is reproduced in
ATSPGM.  A comment card must not contain a block terminator '$'.
Columns 1-5 are used to enter line numbers.  Column 6 is used for
continuation; there is a limit of 19 continuation lines in
FORTRAN.  Columns 7-72, where the block terminator '$' must appear,
contain the statements of the equations.  Columns 73-80 are
ignored.  As in FORTRAN, blanks are not significant.  (A word of
caution, the tab character is not recognized by ATOMCC.)

         The equations in block 1 must be of the form

         DIFF(X,Y,N) = expression,
         variable = DIFF(X,Y,N),
or       variable = expression.

    An expression may contain operations on variables and DIFF(,,)
functions.  The highest order derivative of each dependent variable
must be given explicitly by an equation of the form DIFF(,,) =
expression, but DIFF(,,) functions of lower order may appear in
expressions on the right hand side.  A system of differential
equations may be specified with more than one independent
variable.  But, each independent variable will have the same value
as the solution is computed.  Independent variables are implicitly
defined by the DIFF(,,) function for that independent variable, so
explicitly defining one in an assignment statement will cause an
error message to be printed which indicates that the independent
variable in question has already been defined.

    Incidentally, ATOMCC transforms all constant integer powers of
a variable up to 4 into multiplications.  This is done in order to
avoid the problems raised by the initial value of that variable
being equal to zero.  For integer powers greater than 4, the user
is responsible to see that a l'Hopital situation does not arise.

    Except for integer powers, constant coefficients which appear
in ODE's are converted to real numbers by ATOMCC, so that it does
not matter whether three is written as 3, 3.0, 3.E0, or 3.D0.


                              - 31 -




3.2.2 Parameters in the equations

    Example 3-1 shows a system involving parameters.  It is often
interesting to explore the dependence of the solution on parameters
which appear in the differential equation.  The ATOMCC system is
especially well suited to this problem, because ATSPGM can be
generated only once and then executed repeatedly for different
values of the parameters.

    You may enter the parameters as constants directly into the
statement of the equations.  In that case, ATOMCC writes those
constants directly into ATSPGM; this approach does not allow easy
modification of the parameters.  Note that all constants (except
integer powers) are converted to real numbers by ATOMCC, so that it
does not matter whether three is written as 3, 3.0, 3.E0, or 3.D0.

    If you wish to change the values of the parameters, make each
parameter in the equation to appear as a variable as shown in
Example 3-1. Note that the variable parameters must be assigned
dummy values in block 1, even though the actual values may be
specified at solution time by statements in block 3.

    If the values of the parameters are to be supplied at solution
time by reading from a data file or from a terminal, it may be
convenient to solve the problem repeatedly as shown by Example 3-14
in Section 3.4.2.


3.2.3 COPTION DOUBLE - Double-precision ATSPGM

    Unless you specify differently, the ATSPGM program generated by
ATOMCC is written in single-precision.  A COPTION DOUBLE card, as
the first card in ODEINP, signals ATOMCC to write ATSPGM in
double-precision.  No other changes should be made in block 1. In
particular, if the ODE's have library functions such as SIN, EXP,
etc., ATOMCC automatically inserts DSIN, DEXP, etc.  into ATSPGM.
You should not make those substitutions yourself.

           Example 3-7.  Double precision object program.


COPTION DOUBLE
      DIFF(Y,T,2) = 6*Y*Y + T    $
      $
C Read initial conditions and integration interval
      OPEN(5,FILE='DATA')
      READ(5,1010) START,END,Y(1),Y(2)
 1010 FORMAT(4F10.3)
      WRITE(*,1020) START,END,Y(1),Y(2)
 1020 FORMAT(' Solve the first Painleve transcendent' /
     + ' Interval:          ',2F10.3 /
     + ' Initial conditions:',2F10.3 /)       $
      $


                              - 32 -



In a double-precision ATSPGM program:-

  -  real variables are declared as double precision;

  -  double-precision FORTRAN functions (i.e. DLOG) are generated
     in ATSPGM; (The user must still use single-precision functions
     in the first block.)

  -  call statements in ATSPGM reference the double-precision
     ATOMCC library routines DRSET, DSTEP, and DRDCV;

  -  FORMAT statements use D-format for real variables;

  -  constants are generated in double-precision form.

    Some changes may be required in blocks 2, 3, and 4. ATOMCC
copies them directly into ATSPGM, so you are responsible for any
changes such as declarations or format modifications.

    Be sure to link ATSPGM.OBJ with the double-precision ATOMCC
subroutine library, DRDCV.OBJ.  Incidentally, stiff problems are
solved only in double-precision.


3.2.4 COPTION COMPLX - Complex ATSPGM

    Including a COPTION COMPLX card as the first card in ODEINP
signals ATOMCC to write a single-precision complex ATSPGM.

               Example 3-8.  Complex object program.


COPTION COMPLX
      DIFF(Y,T,2) = 6*Y*Y + T    $
      $
      OPEN(5,FILE='DATA')
      READ(5,1010) MPRINT,NSTEPS,KPTS
 1010 FORMAT(3I5)
      WRITE(*,1010) MPRINT,NSTEPS,KPTS
C
C Read initial conditions
      READ(5,1020) Y(1),Y(2)
 1020 FORMAT(8F10.3)
      WRITE(*,1020) Y(1),Y(2)
C
C Read piecewise linear path
      READ(5,1020) (POINTS(I),I=1,KPTS)
      WRITE(*,1020) (POINTS(I),I=1,KPTS)  $
      $

    No other changes should be made in block 1. If the ODE's have
library functions like SIN, EXP, etc., ATOMCC automatically inserts
CSIN, CEXP, etc. into ATSPGM.  You should not write any complex
functions yourself.

    Some changes may be required in blocks 2, 3, and 4. ATOMCC
copies them directly into ATSPGM, so you are responsible for any


                              - 33 -


changes such as declarations or format modifications.  The
independent and dependent variables are complex, so each must have
both a real and an imaginary part.  So, FORMAT 1010 and 1020 are
changed from Example 3-7 to Example 3-8. We suggest that you study
an example of a complex ATSPGM to see which variables are of type
COMPLEX before you attempt to execute the program.

    The other important change required for a complex ATSPGM is the
manner in which the path of integration is specified.  The path of
integration is a piecewise linear path in the complex plane of the
independent variable.  The path consists of KPTS points stored in
the complex array POINTS as shown in Example 3-8. The complex path
of integration is discussed further in Section 3.4.14.

    For a normal search of the complex plane for singularities, it
is best to keep MPRINT=4.  Higher values of MPRINT only leads to
ever more confusing amounts of print outputs.


3.2.5 COPTION DOUBLE, COMPLX - Double-complex ATSPGM

    Including a COPTION DOUBLE,COMPLX card, as the first card in
ODEINP, signals ATOMCC to write a double-precision complex ATSPGM.

            Example 3-9.  Double complex object program.


COPTION DOUBLE,COMPLX
      DIFF(Y,T,2) = 6*Y*Y + T    $
      $
      OPEN(5,FILE='DATA')
      READ(5,1010) MPRINT,NSTEPS,KPTS
 1010 FORMAT(3I5)
      WRITE(*,1010) MPRINT,NSTEPS,KPTS
C
C Read initial conditions
      READ(5,1020) Y(1),Y(2)
 1020 FORMAT(8F10.3)
      WRITE(*,1020) Y(1),Y(2)
C
C Read piecewise linear path
      READ(5,1020) (POINTS(I),I=1,KPTS)
      WRITE(*,1020) (POINTS(I),I=1,KPTS)  $
      $

    No other changes should be made in block 1. In particular, if
the ODE's have library functions like SIN, EXP, etc., ATOMCC
automatically inserts CDSIN, CDEXP, etc.  into ATSPGM.  You should
not make those substitutions yourself.  Example 3-10 is a portion
of ATSPGM generated for the input shown in Example 3-9.

               Example 3-10.  ATSPGM for Example 3-9


c*+*+*+*+*+*
c     This program was produced by the  ATOMCC  translator version 7.10
c                                       Copyright (C) 1985, Y. F. Chang
c*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+*+


                              - 34 -


c Portions (c) Copyright, Microsoft Corp., 1981.  All rights reserved.
c This program was written for the following inputs
c
COPTION DOUBLE,COMPLX
C     DIFF(Y,T,2) = 6*Y*Y + T
c--------
c no instructions in second input block
c--------
      COMMON /IPASS/ LENSER,LENVAR,MPRINT,MSTIFF,LRUN,
     + KTRDCV,KNTSTP,KTSTIF,KXPNUM,KDIGS,KENDFG,NTERMS,NOPT
     A /RPASS/ RADIUS,ERRLIM,ADJSTF,RCREAL,RCIMAG,RDCERR
     B /CPASS/ START,END,ORDER
     C /DPASS/ H,HNEW
     D /PATHCM/ POINTS,VECTOR,KPTS,KPAST
      COMPLEX*16 START,END,POINTS(10),VECTOR,CZRO,DCMPLX
      COMPLEX ORDER
      DOUBLE PRECISION H,HNEW,AL
      COMPLEX*16 TMPS(  36,  1)
      CHARACTER*6 NAMES
      EQUIVALENCE (TMPS(1,1),Y(1))
      COMPLEX*16 Y(36), T(2), TMPAAB(30), TMPAAA(30)
      COMPLEX*16 SHIFT
      DIMENSION NAMES(1)
      DATA NAMES(1)/'Y.....'/
      Y(33) = 1.1
    9 FORMAT(2X,11HAbove is at,1P2D12.4,9H step no.,I4/)
   10 FORMAT(72H   ATOMCC  Ver. 7.10, Copyright (C) 1985, Y. F. Chang; S
     Aolution results./9H   ******)


      .  .  .  .


      CZRO = DCMPLX(0.D0,0.D0)
c--------
c start of third  input block
c--------
      OPEN(5,FILE='DATA')
      READ(5,1010) MPRINT,NSTEPS,KPTS
 1010 FORMAT(3I5)
      WRITE(*,1010) MPRINT,NSTEPS,KPTS
C
C Read initial conditions
      READ(5,1020) Y(1),Y(2)
 1020 FORMAT(8F10.3)
      WRITE(*,1020) Y(1),Y(2)
C
C Read piecewise linear path
      READ(5,1020) (POINTS(I),I=1,KPTS)
      WRITE(*,1020) (POINTS(I),I=1,KPTS)
c--------
c end of third  input block
c--------


      .  .  .  .


                              - 35 -


c--------
c Calculate radius of convergence and take optimum step.
c--------
         CALL CDRDCV(TMPS,LENVAR,NUMEQS,NAMES)
c--------
c no instructions in fourth input block
c--------
         GO TO (26,28), KENDFG
   26    AL = RADIUS
         H = DSIGN(AL,H)
         IF(MPRINT.GT.4 .OR. NOPT.EQ.0) GO TO 24
         WRITE(*,9) START,KNTSTP
   24    START = START + HNEW*VECTOR
         IF(MPRINT.LT.5) GO TO 27
         WRITE(*,11) KNTSTP,START
         K = Y(33)
         WRITE(*,12) NAMES(K),Y(1), Y(2)
   27 CONTINUE
      WRITE(*,14) NSTEPS
   28 CONTINUE
   29 STOP
      END

    Some changes may be required in blocks 2, 3, and 4. ATOMCC
copies them directly into ATSPGM, so you are responsible for any
changes such as declarations or format modifications.  The
independent and dependent variables are complex, so each initial
condition must have both a real and an imaginary part.

    The important change required for a complex ATSPGM is the
manner in which the path of integration is specified.  The path of
integration is a piecewise linear path in the complex plane of the
independent variable.  The path consists of KPTS points stored in
the complex array POINTS as shown in Example 3-8. The complex path
of integration is discussed further in Section 3.4.14.

    For a normal search of the complex plane for singularities, it
is best to keep MPRINT=4.  Higher values of MPRINT only leads to
ever more confusing amounts of print outputs.


3.2.6 COPTION LENVAR=n - Series length

    The size of the arrays which contain each series is stored in
LENVAR.  The value of LENVAR may be changed by placing the
expression "LENVAR=n" on a COPTION card.  The length of the series
used is controlled by the variable LENSER (see Section 3.4.9).

                  Example 3-12.  COPTION LENVAR=n


COPTION LENVAR=100
      DIFF(POWER,TIME,4) = 3*DP2
      DP2 = DIFF(POWER,TIME,2)        $
      $
      OPEN(5,FILE='DATA')
      READ(5,1010) START,END,(POWER(I),I=1,4)
 1010 FORMAT(6F10.3)


                              - 36 -


      WRITE(*,1010) START,END,(POWER(I),I=1,4)  $
      $

    There are several circumstances under which you may wish to use
a series length which differs from the 30-term series used by
ATSPGM:-

  -  The system of ODE's contains no functions or products of
     dependent variables.  The cost of generating the series is
     then proportional to the series length (instead of the length
     squared), so using longer series may result in a faster
     execution time.

  -  A series begins with many zero terms.  Since ATSPGM requires
     series which have at least 8 nonzero terms, you must lengthen
     the series.

    Section 3.4.9 contains a discussion of the effect of series
length on the execution time of ATSPGM.


3.2.7 COPTION DUMP=n - Diagnostic messages

    This feature involves the operation of the ATOMCC translator
itself and will rarely be used.  It was used during the development
of the software and is documented here only for completeness.  The
form is COPTION DUMP=n, where n indicates the amount of dump to be
written (on your screen).  The following is dumped for each n.

     5   after lexical and syntactical analysis,
and  4   after first optimization,
and  3   after equation sorting and variable type identification,
and  2   after implicit operator and operand analysis,
and  1   after second optimization.

    It is suggested that users who wish to consult with the authors
about a suspected bug in ATOMCC should have available the output
from a run in which "DUMP=5" was specified and directed to a file
for printing.




3.3 Using block 2


    The second, third, and fourth blocks are not processed by
ATOMCC in the same way as the first input block.  ATOMCC merely
copies the data from your ODEINP file directly into ATSPGM.

    The second block is optional and is used to insert
non-executable FORTRAN statements at the beginning of ATSPGM.  This
block gives the user the ability to insert 'SUBROUTINE',
'FUNCTION', 'COMMON', 'DATA', and 'DIMENSION' into ATSPGM.  Example
3-13 shows where statements supplied in block 2 are inserted into
ATSPGM.


                              - 37 -




3.3.1 Subroutine form of ATSPGM

    ATSPGM is normally generated as a main program which can then
be run to solve the specified problem.  The ATOMCC translator has
the capability to generate subroutine or function forms for
ATSPGM.  The subroutine or function can then be called from a user
supplied main program, or from another ATOMCC-generated main
program, to obtain the solution of the problem.

    A subroutine ATSPGM is specified by placing a SUBROUTINE
statement in the second block.  This statement is placed unchanged
at the beginning of the generated code, so you have complete
control of the name of the subroutine and its parameter list.
However, this also means that you have complete responsibility for
passing values to and returning values from the generated program
segment.  ATSPGM written by ATOMCC will have a RETURN statement
instead of a STOP at its end.  Example 3-13 shows the ODEINP for a
subroutine ATSPGM named DIFEQU, which solves the equation f' =
log(sin(x) + f).

              Example 3-13.  Subroutine form of ATSPGM


      DIFF(F,X,1) = ALOG(SIN(X) + F)   $
      SUBROUTINE DIFEQU(COND)        $
      F(1) = COND $
      $

with its calling program

C Driver program for Example 3-13.
C
C A very simple subroutine, START and END are passed through COMMON.
C
      COMMON/CPASS/ START,END,ORDER
 1010 FORMAT(6F10.3)
      READ(5,1010) COND,START,END
      WRITE(*,1010) COND,START,END
      CALL DIFEQU(COND)
      STOP
      END

    Many of the variables used in ATSPGM appear in COMMON, so they
cannot be passed as parameters.  Values can be passed through the
parameter list as temporary variables, or they can be passed in
COMMON by including the appropriate COMMON block in the calling
program as illustrated by Example 3-13.


3.3.2 User declarations

    You may need to declare the type of the additional variables
you introduce into ATSPGM, if you use the double-precision or
complex forms of ATSPGM.  You may include appropriate DIMENSION,
DOUBLE, COMMON, DATA, etc.  in block 2. Any non-executable FORTRAN
statements inserted in the second block are copied into the


                              - 38 -


beginning of ATSPGM exactly as they are written, so you have total
control over the statements entered.  These statements must conform
to FORTRAN specifications; ATOMCC does not check their syntax or
the order of placement.


3.3.3 Common blocks for user

    You may wish to use COMMON blocks.  COMMON declarations may be
included in block 2 like all non-executable FORTRAN statements
discussed in Section 3.3.2.

    Note that many of the variables used in ATSPGM appear in COMMON
blocks and must not be assigned to other blocks.




3.4 Using block 3


    The third block must be used to specify the interval of
integration and the initial conditions.  It may also be used to
change the default values of method-control variables.  Default
values are changed by inserting statements which assign new values
to these variables.  If desired, other statements can be inserted
into ATSPGM at this point.  The user should have an understanding
of the form of ATSPGM and the relative position of the third block
in ATSPGM.  Examples 2-1 and 2-3 and Examples 3-1 and 3-3 are a
pair of ODEINP files with their respective ATSPGM programs.  you
can better understand how the method-control variables work by
studying the statements in the neighborhood of block 3.

    The assignment of values to each of the variables discussed in
this Section may be done by:

  1.  a READ statement (see Examples 2-1, 2-2, 3-1, 3-7, 3-11,
      3-12),

  2.  an assignment statement (see Example 2-6), or

  3.  being passed as temporary variables through a parameter list
      (see Example 3-13), or

  4.  common block shared with a driving main program (see Example
      3-13).


                              - 39 -




3.4.1 Initial conditions

    The initial conditions must be specified in the third block.
The initial values at xo = START of the dependent variable, say
yyy, and its derivatives are assigned to the array yyy as follows:-

         YYY(1) denotes yyy at xo,
         YYY(2) denotes yyy' at xo,
         YYY(3) denotes yyy'' at xo, etc.

    The number of initial conditions which must be specified for
each dependent variable equals the highest derivative of that
variable which appears in the system of equations.  If the system
includes a DIFF(y,x,5), then five initial conditions must be
included for the dependent variable y.  ATOMCC does not check
whether initial conditions have been supplied, since it is possible
that the series variable has been passed through a subroutine
parameter list.  Therefore, it is your responsibility to see that
the required initial conditions are defined.

    When ATSPGM is a subroutine, the series variable may be passed
in the parameter list.  In this case, the initial conditions may be
stored in the proper elements of the series array before ATSPGM is
executed.  An alternative to passing the series via the parameter
list is to pass the initial conditions as temporary variables in
the parameter list.  Assignment statements in the third block then
assign the values to the appropriate dependent variable array
elements.

    Note to those who are familiar with the Taylor series method:-
the initial conditions should be entered as explained above; that
is, Y(3) = y''(xo), not y''(xo)*h*h/2, since this shifting is done
automatically by ATSPGM.


3.4.2 Parameters in the differential equations

    It is often interesting to explore the dependence of the
solution on parameters in the ODE's.  Section 3.2.2 and Example 3-1
showed how the system of equations is entered and how values of the
parameter are assigned.  A refined version of Example 3-1 is shown
in Example 3-14. Here, a loop is used to read successive values of
the parameter.  See the object program listing in Example 3-3 and
observe that statement '28 CONTINUE' is provided by ATOMCC for this
purpose.

             Example 3-14.  Read parameters repeatedly.


      DIFF(X,T,2) = - ALPHA*X*R
      DIFF(Y,T,2) = - ALPHA*Y*R
      R = (X*X + Y*Y)**(-1.5)
      ALPHA = 0.65             $
      $
C Read and echo print code, number of integration steps.
      OPEN(5,FILE='DATA')


                              - 40 -


      READ(5,1010) MPRINT,NSTEPS
 1010 FORMAT(2I10)
      WRITE(*,1010) MPRINT,NSTEPS
C
C Read and echo interval and initial conditions to be used each time.
      READ(5,1020) OLDSTR,OLDEND,X1,X2,Y1,Y2
 1020 FORMAT(6F10.3)
      WRITE(*,1020) OLDSTR,OLDEND,X1,X2,Y1,Y2
C
C Loop for different values of the parameter.
      DO 28 IPROB=1,20
      READ(5,1020) ALPHA
      IF(ALPHA.EQ.0.0) STOP
      WRITE(*,1020) ALPHA
      START = OLDSTR
      END = OLDEND
      X(1) = X1
      X(2) = X2
      Y(1) = Y1
      Y(2) = Y2          $
      WRITE(7,1030) KINTS,START,X(1),X(2),Y(1),Y(2)
 1030 FORMAT(I5,1P5E14.5)   $


3.4.3 Solve a problem repeatedly

    In Example 3-14, the statement 'DO 28 IPROB=1,20' is included
in block 3 to solve the same system of differential equations as
many as 20 times without restarting the program.  The statement '28
CONTINUE' near the end of ATSPGM is written by ATOMCC for this
purpose.  Within this loop, you may vary values of parameters as in
this example, you may vary the initial conditions, or you may vary
the method-control variables.


3.4.4 START, END - Interval of integration

    The integration interval must be specified in the third block
using the two reserved words START and END.  Their use is
illustrated by all the examples.

    If a complex ATSPGM is being used, then the interval of
integration is a piecewise linear path in the complex plane of the
independent variable.  The specification of complex paths of
integration is discussed in Sections 3.4.13 and 3.4.14.


3.4.5 NSTEPS - Number of integration steps, default=40

    The maximum number of integration steps taken during the
solution is the variable NSTEPS.  The default value (40) can be so
small because the use of long series allows very large integration
steps.  You may change this value by inserting a statement in the
third input block which assigns a new value to NSTEPS, or you may
read in a value for NSTEPS from a data file.

    For some stiff problems and in searching for singularities in
the complex plane, it is best to set NSTEPS at a large value.


                              - 41 -




3.4.6 H - Initial trial stepsize, default = 1.0

    The default value of the suggested initial stepsize can be
changed by assigning the desired value to the variable H in the
third block.  The term "suggested initial stepsize" is used because
this value may be adjusted by ATSPGM before the first step is
taken.  This adjustment is made if its need is indicated by an
analysis of the Taylor series for underflow/overflow.  The user can
rarely make a better choice than ATOMCC.

    For stiff problems, in addition to the adjustment of H in
ATSPGM at the first integration step, the stepsize H is adjusted
within DRDCV at every step.  The adjustments at the first integra-
tion step may need some manual assistance.  In such cases, observe
the values generated by the ATSPGM program and project forward to
the next reasonable value and enter it for H.


3.4.7 ERRLIM - Preset accuracy of the solution

    The local error tolerance is the variable ERRLIM.  ATSPGM will
keep the maximum local error less than ERRLIM.  The magnitude of
ERRLIM is automatically set by ATOMCC to be close to the computer
round-off error in both single- and double-precision.  You may set
ERRLIM to be much larger; however, your results will become
inaccurate.  You may not set ERRLIM to be much smaller.

    The error analysis for normal problems is so precise that one
can almost expect the global error to be under control.  This is
one numerical method where the global error is proportional to the
local error.  Therefore, to determine the magnitude of the global
error, one only needs to run the solution a second time with an
ERRLIM set one order of magnitude smaller than the first run and
compare the two results.


3.4.8 ADJSTF - Error control for stiff problems

    The error analysis for stiff problems is not nearly as well
developed as the error analysis for normal problems.  In solving
stiff problems, it is necessary to make assumptions regarding both
the exponential function and the polynomial that are used to fit
the solution being sought.  Therefore, an error-controlling
parameter separate from ERRLIM is used.  The default value for
ADJSTF is 1.E-2. Also, the value of ERRLIM is fixed to 1.E-6 for
stiff problems.  The user is advised to change ADJSTF to lower
values and make additional runs of the solution to be certain of
its correctness.

    There is absolutely no connection between a particular value of
ADJSTF and the magnitude of the error in the computations.  The
most that is claimed is that as the value of ADJSTF is adjusted
smaller, there is likely to be a lowering of the computational
error.  Since the nature of stiff solutions is an approximation,
there can be no true error control.  And so, none is attempted.


                              - 42 -




3.4.9 LENSER - Length of series used, default=30

    The length of the Taylor series is the variable LENSER.  This
value may be changed by assigning a new value to LENSER in the
third block.  However, there are some restrictions governing how
this is done.

    LENSER may be set to any integer between 15 and 30 without any
other changes.  Series of fewer than 15 terms should not be used.
If a series of more than 30 terms is desired, the size of all
series variables (LENVAR) must be increased (see Section 3.2.6).
The user is reminded that the execution time is related to the
length of the series.  The default value of LENSER=30 have been
found to be optmial for fast execution.

    In stiff solutions, LENSER is automatically set to either 15 or
10 depending on the parameter MSTIFF.  In these cases, the user
cannot change these set values of LENSER without going into the
subroutine DRDCV.


3.4.10 MPRINT - Amount of print produced, default=4

    The amount of printout produced by ATSPGM during the solution
of the problem is controlled by the variable MPRINT.  Values of
every dependent variable at each integration step is printed with
an MPRINT=4.  The user can change the default by assigning MPRINT a
different value in the third block.  The amount of printout
produced for values of MPRINT is listed below.

   0  Used for timing purposes; printout is produced only when a
      fatal error occurs.

   1  No print is produced, but the loop controlled by RSET is
      activated to produce user controlled print (see Section 3.5.2).

   2  Print is produced only at points selected by the user.  The
      printout consists of the integration step number, the value of
      the independent variable and the initial conditions for each
      dependent variable.

   4  (default) Print the information for 2 at every integration step.
      (In complex solutions, only the singularity locations are
      printed here.)

   5  In addition to the output under 4, the actual stepsize used
      at each integration step is printed.  In stiff problems, the
      exponential function and the length of the polynomial are
      printed.  In stiff problems, the estimated HSTF and the
      relative goodness of the exponential fit are printed.  (In
      complex solutions, MPRINT=5 is equivalent to 4)

   6  In addition to the information for 5, print (a)the computed
      radius of convergence, (b)location of the singularity(ies),
      and (c)which test was used to locate the poles.


                              - 43 -


   7  In stiff problems, the entire results for the exponential fit
      are printed.

   8  All of MPRINT=6, plus (a)the estimated error in h/Rc, (b)the
      values of the elements of each series, and (c)messages for
      all tests used for the radius of convergence.

  10  All diagnostic messages including all the series terms.

    Users who will use ATOMCC often should try setting MPRINT=10 to
become familiar with the information available during the solution
of a problem.

    The value of MPRINT may be dynamically controlled during the
computation of a solution by inserting statements in the fourth
block which test the current value of START (or KINTS), and set the
value of MPRINT accordingly.


3.4.11 DLTXPT - Print point increments, default = 0.0

    ATSPGM uses variables XPRINT and DLTXPT to control values of
the independent variable for which print is produced.  We will
discuss a variety of ways these variables can be used.

    If DLTXPT=0.0, then print is produced only at the integration
steps chosen by the program.  This is the default condition.

    Print is produced at equally spaced points by specifying DLTXPT
= xxx (the desired spacing) as shown in Example 2-8 in Section 2.3.
A statement may be placed in the third block which assigns the
desired value to DLTXPT.  You should also specify MPRINT=2.
Otherwise print is produced both at your selected points and also
at the integration steps selected by ATSPGM.  MPRINT is discussed
in Section 3.4.10.

    The values of each dependent variable are printed at the
equally spaced points by expanding a series which has already been
computed.  The integration does not step to the equally spaced
points.  Hence, requesting intermediate output has no effect on the
number or size of integration steps taken.

    More creative print points are possible by the use of block 4.
See Example 3-15 in Section 3.5.3.

    DLTXPT can be used in stiff solutions to generate output at
desired values of the independent variable.  DLTXPT cannot be used
in conjuction with ZEROT, where the solution may be stopped for any
value of any function.


3.4.12 KTRDCV - Dynamic suppression of CALL RDCV

    KTRDCV can be used to speed up the execution of ATSPGM for
systems of ODE's for which some components of the solution do not
constrain the integration stepsize.  For KTRDCV=N, the series
analysis is performed only for the first N components of the
solution.  You can generate ATSPGM once, run it to see the radii of


                              - 44 -


convergence of each component, and use KTRDCV on subsequent runs to
reduce the execution time of ATSPGM.  You should order the input
ODE's such that those with larger radii are listed last.  Careful
study of the ATSPGM program for a sample system of ODE's will help
you see how KTRDCV works.


3.4.13 KPTS - Number of points on complex path

    In Section 3.4.4, we discussed the specification of the
interval of integration using START and END.  If ATOMCC has been
directed to generate complex object code (COPTION COMPLX, Sections
3.2.4 or 3.4.5), then the path of integration is a piecewise linear
path in the complex plane of the independent variable.  See Example
3-8 in Section 3.2.4.

    The variable KPTS is the number of vertices belonging to that
piecewise linear path, including both of the endpoints.  The
complex array 'POINTS' holds the vertices.  POINTS(1) becomes
START, and POINTS(KPTS) becomes END.  The value of the solution is
printed at each element of POINTS.  KPTS may be at most 10.


3.4.14 POINTS - Complex path of integration

    The complex array 'POINTS' specifies the path of integration in
the complex plane of the independent variable.  Its use with KPTS
is discussed in Section 3.4.13. There may be at most 10 points
specified.


3.4.15 MSTIFF=10 - Solutions which are entire

    Solutions which are entire (have no singularities in the finite
plane), should not be solved using the ATOMCC system.  It is a
total waste of computing power to solve linear problems using
ATOMCC.  This is particularly true for linear 'stiff' problems.

    It can be EASILY show that ALL solutions that are entire can be
solved in quasi-closed forms.  This INCLUDES two-point boundary
value problems!

    If ATOMCC recognizes that a system of ODE's involves no
functions or products (an entire solution), it sets MSTIFF=10 in
ATSPGM.  Subroutine RDCV recognizes MSTIFF=10 as a flag for a
special test for series that is entire.


3.4.16 MSTIFF=20,21,22 - Stiff problems.

    This version of ATOMCC contains a double-precision algorithm to
solve stiff problems.  To use it, one can either set MSTIFF=20, or
21, or 22. Other parameters that should be controlled are H,
ADJSTF, and NSTEPS.  It is also desirable to set MPRINT to 7, at
least initially, for observing the progress of the solution.  If it
should be evident that the problem is not really stiff, then it is
most advisable to solve it as a normal problem.


                              - 45 -


    MSTIFF=20 is the more conservative of the three algorithms.  In
this case, LENSER is set to be 15. The default value for ADJSTF,
the error-controlling parameter, is a rather large 1.E-2. The user
should run the stiff solution at least one more time with a
somewhat smaller ADJSTF, say 1.E-3, to check on its validity.
Since the default value for NSTEPS is 40, this should definitely be
set at 500 or more, to be sure that the stiff solution can be
completed.  The initial stepsize H may need some adjustment by the
user.  He should study the H-adjustment messages from ATSPGM and
take over control if and when the automatic adjustment is incapable
of reaching a desirable value for H. When MSTIFF=20, those stiff
problems that have steady-state solutions are identified.  After
which, one should perform some manual manipulations before
re-submitting the problem to ATOMCC.

    When MSTIFF=21, LENSER is set to only 10. So, this option
should be used only if the user is absolutely certain that the
problem under study is very stiff.  The solution of stiff problems
under this option is considerably faster than that for MSTIFF=20,
because not only are the series shorter, the integration stepsizes
are considerably larger.  The same statements as above applies for
the parameters H, ADJSTF, and NSTEPS.  (One of the Enright-Hull set
of stiff problems, E2, was solved by ATOMCC in one step!)  There is
no attempt to identify steady-state solutions when MSTIFF=21.

    MSTIFF=22 is identical to MSTIFF=20 except for the fact that
there is no attempt to identify steady-state solutions.

    As mentioned above, the stiff algorithm is written in double-
precision.  It is simply not cost effective to solve such problems
using single-precision.  There is one other restriction on stiff
problems.  All such problems must be stated as first-degree ODE's.

    The regular printing option, DLTXPT, functions properly in
stiff solutions.  So, it is possible to obtain uniformly spaced
print points or logarithmicaly spaced print points.


3.4.16.1 Steady-State Stiff Problems

    There are some stiff problems that approach steady-state solu-
tions.  Sometimes this occurs with all of the functions becoming
constant simultaneously.  In such a case, the ATOMCC solution will
stop and points out the fact that every function seem to have
reached constant values.  The user can decrease ADJSTF and repeat
the solution to verify the truth of this fact.  He can also run the
problem with MSTIFF=22 and observe the perpetual singleness of the
results.

    In other instances, a particular function reaches constancy
before any of the others.  When this happens, the following message
will be printed, with obvious meaning.

The function Y4.... is constant at  3.115283D-04
  Look for a steady-state solution.
  Set all derivatives to zero, use the value of Y4.... given above,
    and solve for the other functions in easy situations.
  Then, re-submit to ATOMCC.  Use MSTIFF=21 or 22, do not use MSTIFF=20.


                              - 46 -



    This example is for the problem CHEM6 out of the Enright-Hull
set of chemical stiff problems.  This solution stopped at time =
0.0561, where the solution had hardly gotten started.  After
solving the easy cases as suggested, the ATOMCC solution of the
re-submitted problem reached the final time of 1000 in 26 steps
using MSTIFF=21!

    The equational input for the first ATOMCC run is as follows.

      DIFF(Y1,T,1) = 1.3*(Y3-Y1) + 10400*AK*Y2
      DIFF(Y2,T,1) = 1880*(Y4 - Y2*(1+AK))
      DIFF(Y3,T,1) = 1752 - 269*Y3 + 267*Y1
      DIFF(Y4,T,1) = 0.1 + 320*Y2 - 321*Y4
      AK = EXP(20.7 - 1500/Y1)

The equational input for the second ATOMCC run, after encountering
the constancy message is as follows.

      DIFF(Y1,T,1) = 1.3*(Y3-Y1) + 10400*AK*Y2
      DIFF(Y2,T,1) = 1880*(Y4 - Y2*(1+AK))
      Y3 = 1752/269 + 267/269*Y1
      Y4 = 3.115283E-04
      AK = EXP(20.7 - 1500/Y1)

The differences between to two inputs are in the third and fourth
equations.




3.5 Using block 4


    ATOMCC copies the fourth block directly into ATSPGM near the
end of the code for each integration step.  Example 3-3 in Section
3.1.3 shows the location of block 4 in ATSPGM.  The fourth block is
used to tailor ATSPGM to your special requirements.  Several
examples are provided in Section 3.1.1, but many other creative
uses are possible.  Most applications require a good knowledge of
the ATOMCC system for solving ODE's.


3.5.1 Automatic printing of output points

    The ATSPGM program automatically prints the values of each
component of the solution at each integration step (see Example 2-7
in Section 2.2.6). It can print the same information at equally
spaced points you select using DLTXPT (see Section 3.4.11). This
technique for generating output at equally spaced points requires
no use of block 4, but it will help you to understand subsequent
sections if you understand how ATSPGM generates this equally spaced
output.  Please refer to the object code in Example 3-3 in Section
3.1.3 as you read.

    The points selected for output do not affect the integration
steps so that the local error remains proportional to the global
error.  For each output point within an integration step, the


                              - 47 -


solution is computed by evaluating its Taylor series.  Inside the
subroutine RSET, each output point within the the integration step
is controlled by a loop which begins with the statement '24 IF
(KENDFG.EQ.3) KENDFG=1' and ends with the statement 'GO TO
(26,28,24), KENDFG'.  KENDFG is a flag which controls the loop.

    Subroutine RDCV sets KENDFG=2 if the current steps reaches END,
otherwise KENDFG=1.  Subroutine RSET determines whether the next
output point (XPRINT) lies within the next step of integration.  If
not, RSET leaves KENDFG unchanged (1 or 2) and stores the initial
conditions required by the next step.  If there is a print point
within the next integration step, RSET prints the solution at that
output point (suppressed if MPRINT=0 or 1), stores values of the
series and its derivatives as elements LENSER+1, LENSER+2,...  in
the series, and sets KENDFG=3.  Then, the 'GO TO' in ATSPGM returns
to label 24 to handle additional output points within the next
integration step.  Once all of the output points within the next
step have been passed, RSET returns KENDFG to its original value (1
or 2), and the solution proceeds forward.

    You can use DLTXPT to generate output at equally spaced points
without understanding how that output is generated, but if you wish
to produce some special output as discussed in the following
Sections, this understanding is necessary.


3.5.2 User controlled printing of output points

    The object program automatically prints the values of each
component of the solution at each integration step (see Example 2-7
in Section 2.2.6). The method used to produce output at user
selected points is discussed in Section 3.5.2. However, you may use
block 4 as shown in Example 3-1 in Section 3.1.1 to produce output
according to your particular needs.

    The execution time for ATSPGM is sensitive to the amount of
output produced.  Unless you are interested in the automatically
produced output, place MPRINT=1 in block 3. Then, RSET controls the
necessary looping as discussed in Section 3.5.1, but it produces no
output.  This yields the solution at equally spaced points.
Unequal spacing can be achieved by changing DLTXPT within block 4.
Example 3-15 in Section 3.5.3 shows how to obtain the solution at
logarithmically spaced points.

    Many variables are accessible for your use in block 4. The
elements of the array for each dependent variable, [y(LENSER+1),
y(LENSER+2), etc.], contain the values of that variable and its
derivatives evaluated at XPRINT.  The user can therefore print
these values totally independent from the print produced by
ATSPGM.  Additional usable values are stored in the top three
positions of the series; y(LENVAR-2) is the series length actually
used for that particular variable, y(LENVAR-1) is the multiplying
constant for the exponential function in stiff solutions, and
y(LENVAR) is the exponential coefficient in stiff solutions.


                              - 48 -




3.5.3 Logarithmic spacing of output points

    The techniques discussed in Sections 3.5.1 and 3.5.2 produce
values of the solution at equally spaced points.  Output at
non-uniformly spaced points is obtained by adjusting DLTXPT and
XPRT3 within block 4 as shown by Example 3-15. The print points are
at r = 100, 50, 20, 10, 5, 2, etc.  This example also shows that
problems can be integrated in the negative direction.

                 Example 3-15.  Logarithmic print.


      DIFF(F,R,2) = (F**3 - F)/R**2  $
      $
      START = 100.0
      END = 1.0E-8
      F(1) = 0.99
      F(2) = 1.0E-4
      DLTXPT = -50.0
      NPTR = 1
      $
      IF(KENDFG .NE. 3) GO TO 25
      MPRINT = 2
      GO TO (501,503,502), NPTR
  501 DLTXPT = - 0.6*XPRINT
      GO TO 504
  502 NPTR = 0
  503 DLTXPT = - 0.5*XPRINT
  504 NPTR = NPTR + 1
      XPRT3 = XPRINT + DLTXPT  $


3.5.4 ZEROT - Stopping and printing at roots of variables

    It is often of interest to locate points at which a component
of the solution has a root or assumes some specified value.  The
subroutine ZEROT, in the ATOMCC subroutine library in single-
precision, do automatically solve such problems.  DZEROT is the
double-precision version.

    The form of the CALL is

             CALL ZEROT(NUMBER,Y,ROOT,KEY,TMPS,LENVAR,NUMEQS)

    where

    NUMBER is the index of the Y series term whose root is sought,
    Y      is the variable whose root is desired,
    ROOT   is the value Y is to assume (= 0 for a root),
    KEY    is 1 if Y is a dependent variable, or
              0 if Y is not a dependent variable.

    The arguements TMPS, LENVAR, and NUMEQS must be exactly as
written above.


                              - 49 -



               Example 3-16.  Rootfinding with ZEROT.


      DIFF(Y,T,2) = 6*Y*Y + T  $
      $
      START = 0.0
      END = 1.15
      ROOT = 20.0
      Y(1) = 1.0
      Y(2) = 0.0
      $
      CALL ZEROT(1,Y,ROOT,1)
      $

    When the variable whose root is being sought is not a dependent
variable, the following statements are used.  Note that KEY is set
to 0.

        CALL ZEROT(2,VARY,0.0,0)
        IF(LRUN.NE.0) GO TO 25
        TEMP = START + HNEW
        WRITE(7,1010) KINTS,TEMP,VARY(1),VARY(2)
   1010 FORMAT(I5,3F10.4)
        GO TO 25  $

    In these examples, it is not necessary for one to print the
information as shown in the second case.  The ATSPGM program does
stop and restart the solution automatically at the exact root and
the output is controlled by MPRINT.

    ZEROT requires that the series for the variable be convergent
over the range of the integration stepsize being used.  If the
series is not convergent over this range, then ZEROT cannot
possibly locate the root with any degree of accuracy.  There are
two instances where this convergence requirement is not met.  One
is where the printing steps have been placed under the control of
DLTXPT, and the other is where a stiff problem is being solved.
Therefore, ZEROT cannot be used for a non-zero DLTXPT, and in stiff
problems.

    The index NUMBER can be set at any positive (non-zero) integer
value; however, obviously when NUMBER is very large the accuracy of
the root will suffer.


3.5.5 Finding singularities in real solutions

    (This information is only for those users who are interested in
finding the conjugate pairs of singularities closest to the real
axis.  For a more detailed study of singularities, the user should
invoke COPTION COMPLX, Section 3.2.4)

    A unique feature of the ATOMCC system is its ability to provide
analytic information about the solution.  For example, subroutine
RDCV estimates the location and order of primary singularities at
each integration step in order to compute the optimal stepsize.
This information may be printed using MPRINT=6 (see Section


                              - 50 -


3.4.10). However, several steps are usually required to pass
between a conjugate pair of singularities, so several different
estimates for each location are written.  Estimates made close to
the singularities are more accurate than estimates made from
further away.

    Example 3-17 prints only one estimate for the location and
order of each conjugate pair of singularities.  That estimate is
written on the step which has just begun to recede from the
estimated location.  Hence that step is one of the two steps
closest to the singularity pair.

       Example 3-17.  Finding singularities in real solutions


COPTION DOUBLE
C
C Test dynamic printing of singularity positions.
C
      DIFF(Y1,T,1) = 2*(Y1 - Y1*Y2)
      DIFF(Y2,T,1) = Y1*Y2 - Y2     $
      $
      ERRLIM = 1.0E-6
      START = 0.0
      END = 20.0
      Y1(1) = 1.0
      Y2(1) = 3.0
      MPRINT = 10
      PRTSNG = START
      SNGTOL = 0.1
      WRITE(8,2010)
 2010 FORMAT(///6H  Step,9X,1Hx,10X,2HRc,8X,4HReal,
     A 8X,4HImag,7X,5Horder,7X,5Herror/)
      $
C
C Print locations of singularities.
C
      WRITE(*,2010)
      WRITE(*,2020) KINTS,START,RADIUS,RCREAL,RCIMAG,ORDER,RDCERR
C
C If this is a user print step, then continue.
      IF(KENDFG.EQ.3) GO TO 25
      IF(KINTS.EQ.NSTEPS) GO t0 102
C If this is a normal step, then jump.
      IF(KENDFG.EQ.1) GO TO 105
C Handle the last step  -  we might be approaching a singularity
C  if we have already printed this primary singularity, then continue.
  102 IF(ABS(RCREAL-PRTSNG).LT.SNGTOL) GO TO 25
C If RDCV was uncertain, then continue.
      IF(RDCERR.GE.1.0) GO TO 25
      WRITE(8,2020) KINTS,START,RADIUS,RCREAL,RCIMAG,ORDER,RDCERR
 2020 FORMAT(I4,1P6E12.3)
      GO TO 25
C
C Handle a normal step.
C
C If confused, then continue.
  105 IF(RDCERR.GE.1.0) GO TO 25


                              - 51 -


C If approaching the primary singularity, then continue.
      IF((START-RCREAL)*H.LE.0.0) GO TO 25
C If already printed, then continue.
      IF(ABS(RCREAL-PRTSNG).LT.SNGTOL) GO TO 25
      WRITE(8,2020) KINTS,START,RADIUS,RCREAL,RCIMAG,ORDER,RDCERR
      PRTSNG = RCREAL   $

                  Example 3-18.  Printed locations


Step     x         Rc       Real       Imag      order      error

 1  0.000E+00  5.978E-01 -3.364E-01  4.941E-01  9.457E-01  5.588E-04
 9  5.180E+00  4.954E-01  5.151E+00  4.945E-01  9.784E-01  7.081E-04
19  1.083E+01  5.297E-01  1.063E+01  4.943E-01  9.595E-01  2.888E-04
28  1.627E+01  5.145E-01  1.612E+01  4.943E-01  9.662E-01  2.660E-04


3.5.6 Stopping short of a singularity

    If one component of the solution has a singularity inside the
interval of integration, then the problem does not have a solution
on that interval, so the integration process should stop.  ATSPGM
stops when NSTEPS steps have been taken (see Section 3.4.5), or it
stops when the integration stepsize is so small relative to the
current point of expansion that the solution would not advance.




3.6 Editing of ATSPGM


    While the use of blocks 2, 3, and 4 to insert FORTRAN code
directly into ATSPGM gives you a very powerful and flexible tool,
you may have needs which can only be met by editing ATSPGM
yourself.  Section 3.6.1 uses the relatively common problem of
producing efficient output at specified points to illustrate the
approach.


3.6.1 TERM - Fast generation of printing at output points

    The printing of solution values at equally spaced points has
already been discussed in Sections 2.3, 3.4.11, 3.5.1, and 3.5.2.
The technique discussed in Section 3.5.2 can be used to generate
nearly all of the information about the solution at any points you
select.  The execution time is much faster using MPRINT=1, to
suppress printing by RSET, than with MPRINT=2 or 4.

    The execution time of ATSPGM can be further reduced by
controlling the looping for each print point within the current
integration step by calling the subroutine TERM to evaluate the
series for the solution at any desired point.  This subroutine is
listed as Example 3-19.


                              - 52 -



                 Example 3-19.  TERM source listing


C  SUBROUTINE TERM(SERIES,NTERMS,T,WORKAR)
C
C  ATOMCC-LIBRARY, COPYRIGHT  (C)  1983.
C
C  EVALUATES SERIES AND ITS DERIVATIVES AT T
C
C  AT ENTRY, 'SERIES' CONTAINS A SERIES EXPANDED AT SOME POINT
C  'START' WITH A STEP OF 'H'.
C
C  AT EXIT, WORKAR(1)  = FUNCTION EVALUATED AT T
C         WORKAR(2)  = DERIVATIVE OF FUNCTION EVALUATED AT T
C         . . .
C         WORKAR(NTERMS) = NTERMS - 1 DERIVATIVE EVALUATED AT T
C
C  PARAMETERS:
C    SERIES  -  SERIES BEING SUMMED
C    NTERMS  -  NUMBER OF INITIAL CONDITIONS NEEDED
C    T      -  POINT AT WHICH SERIES IS TO BE EVALUATED
C    WORKAR  -  WORK AREA, SAME TYPE AS SERIES, DIMENSION NTERMS
C
     SUBROUTINE TERM (SERIES, NTERMS, T, WORKAR)
C********  DATA DECLARATIONS
     DIMENSION SERIES(1), WORKAR(1)
     COMMON /CPASS/START,END,ORDER
    A  /DPASS/H,HNEW,XPRINT,DLTXPT
    +  /IPASS/LENSER,LENVAR,MPRINT,MSTIFF,LRUN,KTRDCV,
    +   INTSTP,KTSTIF,KXPNUM,KDIGS,KENDFG,NTERMS
C********  COMPUTE INITIAL CONDITIONS AND STORE THEM IN WORK AREA
     RATIO = (T-START)/H
     DUMMY = H/FLOAT(LENSER)
     DO 40 I=1,NTERMS
      ILEN = LENSER - I - 1
      DUMMY = DUMMY*FLOAT(LENSER - I + 1)/H
      DLOOP = DUMMY
      SUM = SERIES(LENSER)*RATIO*DLOOP
      DO 30 J=1,ILEN
       LMJ = LENSER - J
       DLOOP = DLOOP*FLOAT(LMJ-I+1)/FLOAT(LMJ)
       SUM = (SUM + DLOOP*SERIES(LMJ))*RATIO
   30  CONTINUE
      WORKAR(I) = SUM + DLOOP*SERIES(I)/FLOAT(I)
   40 CONTINUE
C********  RETURN
 9999 RETURN
     END

    To illustrate the use of subroutine TERM as an example of
editing ATSPGM, Example 3-20 shows the ODEINP file which was used
to generate ATSPGM listed as Example 3-21. Three arrays have been
dimensioned to hold the results computed by TERM, and new variables
DXP and FIRSTX have been defined to fill roles analogous to DLTXPT
and XPRINT.


                              - 53 -



               Example 3-20.  ODEINP for TERM example


     S = 10.0
     DIFF(X,T,3) = Y - E*S*X
     DIFF(Y,T,3) = - X*Z + X - E*Y
     DIFF(Z,T,3) = X*Y - E*B*Z
     E = 0.02
     B = 8.0/3.0  $
     DIMENSION USER1(3),USER2(3),USER3(3)  $
     OPEN(7,FILE='DATA')
     OPEN(8,FILE='FAST',STATUS='NEW')
     READ(7,*) X(1),Y(1),Z(1)
     READ(7,*) X(2),Y(2),Z(2)
     READ(7,*) X(3),Y(3),Z(3)
     READ(7,*) START,END,DLTXPT
     READ(7,*) DXP
     FIRSTX = START
     WRITE(*,1010) X(1),Y(1),Z(1),START,END,DLTXPT,FIRSTX,END,DXP
 1010 FORMAT(3F10.4)
     NSTEPS = 5
     WRITE(8,1020)
 1020 FORMAT(3X,1Ht,5X,1Hx,7X,2Hx',6X,3Hx'',
    + 5X,1Hy,7X,2Hy',6X,3Hy'',5X,1Hz,7X,2Hz',6X,3Hz''/)    $
     $

with data file

  1.0,3.0,5.0
  0.0,0.0,0.0
  -1.0,0.0,1.0
  0.0,5.0,0.5
  0.5

    The following lines of code must be inserted just below the
CALL RDCV statement inside ATSPGM.  A portion of the results are
shown in Example 3-22.

           Example 3-21.  Object program for Example 3-20


         CALL RDCV(TMPS,LENVAR,NUMEQS)
c=========== User must insert these lines by editing.
  100 IF(FIRSTX.GE.START+HNEW) GO TO 24
      CALL TERM(X,3,FIRSTX,USER1)
      CALL TERM(Y,3,FIRSTX,USER2)
      CALL TERM(Z,3,FIRSTX,USER3)
      WRITE(8,1030) FIRSTX,USER1,USER2,USER3
 1030 FORMAT(1X,F3.1,1P9E8.1)
      FIRSTX = FIRSTX + DXP
      GO TO 100
c=========== end of user inserted lines.
   24    IF(KENDFG.EQ.3) KENDFG = 1


                              - 54 -



               Example 3-22.  Output for Example 3-20


  t    x       x'      x''     y       y'      y''     z       z'     z''

0.0 1.0E+00 0.0E+00-1.0E+00 3.0E+00 0.0E+00 0.0E+00 5.0E+00 0.0E+00 1.0E+0
0.5 9.3E-01-1.5E-01 3.9E-01 2.9E+00-5.0E-01-2.0E+00 5.1E+00 8.3E-01 2.3E+0
1.0 9.6E-01 3.7E-01 1.6E+00 2.3E+00-2.0E+00-4.1E+00 5.9E+00 2.3E+00 3.4E+0
1.5 1.4E+00 1.4E+00 2.3E+00 6.9E-01-4.8E+00-7.4E+00 7.5E+00 4.2E+00 4.1E+0




3.7 Large systems


    As supplied to you, the ATOMCC translator can handle up to 900
equations.  If you should have a need to increase this limit, a
special program can be easily prepared.




3.8 Solving ODE's in the complex domain


    The ATOMCC compiler allows for the solution of ODE's in the
complex domain.  This unique capability can be used to explore the
structure of the singularities in the complex domain of non-linear
problems.  Linear problems, of course, have entire solution
functions and therefore do not have any singularities of interest
in the finite complex plane.  Non-linear problems may have
singularities which cover the entire complex plane.

    There are essentially two types of non-linear problems; those
with definite limit cycles, and those with strange attractors.  For
the former, the singularities form a regular lattice in the complex
plane.  For the latter, the singularities form structures that defy
simple descriptions.  The purpose of solving ODE's in the complex
domain is to study the structures formed by the singularities.  The
ATOMCC compiler is well suited for this task, and it is the only
method extant that is capable of calculating the precise location
and order of all the singularities in a finite region of the
complex plane of an ODE solution.

    It is simple to use the ATOMCC compiler to search for the
singularities.  First, you must insert a COPTION COMPLX card as the
first card in ODEINP.  This will cause the ATOMCC compiler to
generate ATSPGM that will solve the ODE using paths into the
complex domain.  Secondly, you must specify the path to be taken by
the solution.  This path is fixed by specifying the vertices of
straight-line segments in the path.  The path taken must be
composed of straight-line segments.  The first vertex is of course
the starting point of the solution.  A maximum of ten vertices may
be specified.  These vertices are to be placed into an array called
POINTS, and the number of vertices used is stored in the variable


                              - 55 -


KPTS.  The ATOMCC solution will follow the path thus specified and
locate all the singularities near this path.

    Since the ATOMCC solution will find all the singularities near
the path specified by the user, certain problems may occur.  First,
the user may have by accident specified a path exactly midway
between two singularities.  In this event, there will be no
information output from ATOMCC.  The path must be slightly closer
to one singularity than another; otherwise, ATOMCC cannot find the
nearest singularity.  Secondly, the user may have by accident
specified a path that is too close to a singularity, or perhaps
even a path that is directly pointed at a singularity.  In this
event, the ATOMCC solution will grind away and take very small
steps.  The information from ATOMCC beyond any such close encounter
will be unreliable.  In all cases, the user is well advised to
change the path in the complex plane ever so slightly and make a
second run to double check his results.  In our experience, it is
best to perform the complex integrations using double precision.
Insert a COPTION DOUBLE,COMPLX card as the first card in your
input.  With just minimal experience, the user will find the use of
ATOMCC in the complex plane to be fast and easy.

    For good results, the first leg of the path into the complex
domain should be directed straight up in the imaginary direction.
Do not make the first leg of the path coincident with the real
axis.  This introduces subtraction errors into the complex
solution.

    When there are complex constants in the equations, the user is
entirely responsible for seeing to it that those constants are
properly specified with TYPE declarations in block 2, and the
values are properly entered as CMPLX(--,--) or DCMPLX(--,--) in
block 3. The reason for this requirement on the part of the user is
because if the facts are otherwise, the user may then be required
to use an editor to delete some possibly incorrect specification
written by ATOMCC.  It is better to be a bit short and correct than
to be long and in error.